In this script, we are building the spatiotemporal prediction model for black carbon in the Denver metro area. This model uses data from 5 sampling campaigns (as detailed in 15_Exploring_BC_Data.R). Two models are developed here: one which uses one temporal basis function (Model A) and one that uses two temporal basis functions (Model B).
Updates to this model are in response to reviewer comments on an earlier manuscript version submitted to ES&T.
A summary table at the end of the document summarizes the performance for each of these models. Overall, Model B performs well and allows us to predict BC concentrations at a weekly temporal resolution from 2009 through 2020.
## Loading required package: Matrix
Here I’m reading in the filter data sets and getting them into the right shape for the SpatioTemporal package. Here the ST covariates cover the entire study period (2009 - 2020).
lur_data_all <- read_csv(here::here("Data", "Final_BC_Data_Set_for_ST_Model.csv"))
## Parsed with column specification:
## cols(
## .default = col_double(),
## site_id = col_character(),
## filter_id = col_character(),
## campaign = col_character(),
## StartDateTimeLocal = col_logical(),
## EndDateTimeLocal = col_logical(),
## sample_week = col_date(format = ""),
## logged_runtime = col_logical(),
## logged_rt_volume_L = col_logical(),
## low_volume_flag = col_logical(),
## ultralow_volume_flag = col_logical(),
## pm_below_lod = col_logical(),
## pm_below_loq = col_logical(),
## bc_below_lod = col_logical(),
## negative_pm_mass = col_logical(),
## potential_contamination = col_logical(),
## pm_mass_ug = col_logical(),
## bc_mass_ug_corrected = col_logical(),
## blank_corrected_bc = col_logical(),
## bc_blank_mean = col_logical(),
## pm_ug_m3_raw = col_logical()
## # ... with 36 more columns
## )
## See spec(...) for full column specifications.
lur_data <- lur_data_all
bc_obs <- lur_data %>%
dplyr::select(site_id, st_week, bc_ug_m3) %>%
group_by(site_id, st_week) %>%
summarize(bc_ug_m3 = mean(bc_ug_m3, na.rm = T)) %>%
mutate(bc_ug_m3 = ifelse(is.nan(bc_ug_m3), NA, bc_ug_m3)) %>%
mutate(log_bc = log(bc_ug_m3)) %>%
dplyr::select(-bc_ug_m3)
#' Pivot to a wide data frame
bc_obs2 <- bc_obs %>%
pivot_wider(id_cols = st_week, names_from = site_id, values_from = log_bc) %>%
# names_from = site_id, values_from = bc_ug_m3) %>%
as.data.frame() %>%
arrange(st_week)
rownames(bc_obs2) <- bc_obs2$st_week
bc_obs2$st_week <- NULL
bc_obs2 <- as.matrix(bc_obs2)
bc_weeks <- rownames(bc_obs2)
tail(bc_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
#' Spatial covariates (scaled)
load(here::here("Results", "BC_LASSO_Model_5C.Rdata"))
lasso_coef_df <- data.frame(name = log_bc_lasso_coef3@Dimnames[[1]][log_bc_lasso_coef3@i + 1],
coefficient = log_bc_lasso_coef3@x)
covars <- as.character(lasso_coef_df$name[-1])
covars
## [1] "impervious_2500" "open_2500" "med_int_2500" "high_int_100"
## [5] "dist_m_cafos" "aadt_100" "aadt_250"
covar_fun <- paste("~", paste(covars, collapse = " + "))
covar_fun
## [1] "~ impervious_2500 + open_2500 + med_int_2500 + high_int_100 + dist_m_cafos + aadt_100 + aadt_250"
bc_sp_cov <- select(lur_data, site_id, lon, lat, all_of(covars)) %>%
distinct() %>%
mutate_at(.vars = vars(all_of(covars)),
scale)
summary(bc_sp_cov)
## site_id lon lat impervious_2500.V1
## Length:69 Min. :-105.2 Min. :39.57 Min. :-1.9791101
## Class :character 1st Qu.:-105.0 1st Qu.:39.72 1st Qu.:-0.7093639
## Mode :character Median :-104.9 Median :39.75 Median : 0.1278296
## Mean :-104.9 Mean :39.76 Mean : 0.0000000
## 3rd Qu.:-104.8 3rd Qu.:39.79 3rd Qu.: 0.4803932
## Max. :-104.7 Max. :39.95 Max. : 1.8291463
## open_2500.V1 med_int_2500.V1 high_int_100.V1
## Min. :-1.3794337 Min. :-2.1968789 Min. :-0.532437
## 1st Qu.:-0.5497931 1st Qu.:-0.7166467 1st Qu.:-0.532437
## Median :-0.1781892 Median : 0.2108832 Median :-0.532437
## Mean : 0.0000000 Mean : 0.0000000 Mean : 0.000000
## 3rd Qu.: 0.4072220 3rd Qu.: 0.7753690 3rd Qu.: 0.015894
## Max. : 2.7693377 Max. : 1.7719166 Max. : 4.768092
## dist_m_cafos.V1 aadt_100.V1 aadt_250.V1
## Min. :-2.1610758 Min. :-0.390807 Min. :-0.664113
## 1st Qu.:-0.4079546 1st Qu.:-0.390807 1st Qu.:-0.664113
## Median :-0.0010458 Median :-0.390807 Median :-0.365247
## Mean : 0.0000000 Mean : 0.000000 Mean : 0.000000
## 3rd Qu.: 0.3418171 3rd Qu.:-0.115098 3rd Qu.: 0.063484
## Max. : 2.2925064 Max. : 4.632471 Max. : 4.155107
bc_sp_cov <- bc_sp_cov %>%
rename(ID = site_id) %>%
mutate(lon2 = lon, lat2 = lat) %>%
st_as_sf(coords = c("lon2", "lat2"), crs = ll_wgs84) %>%
st_transform(crs = albers)
sp_coords <- do.call(rbind, st_geometry(bc_sp_cov)) %>% as_tibble()
names(sp_coords) <- c("x", "y")
bc_sp_cov <- bind_cols(bc_sp_cov, sp_coords) %>%
st_set_geometry(NULL) %>%
as.data.frame()
head(bc_sp_cov)
bc_sp_cov$type <- ifelse(bc_sp_cov$ID == "central", "central", "dist")
bc_sp_cov$type <- as.factor(bc_sp_cov$type)
summary(bc_sp_cov)
## ID lon lat impervious_2500.V1
## Length:69 Min. :-105.2 Min. :39.57 Min. :-1.9791101
## Class :character 1st Qu.:-105.0 1st Qu.:39.72 1st Qu.:-0.7093639
## Mode :character Median :-104.9 Median :39.75 Median : 0.1278296
## Mean :-104.9 Mean :39.76 Mean : 0.0000000
## 3rd Qu.:-104.8 3rd Qu.:39.79 3rd Qu.: 0.4803932
## Max. :-104.7 Max. :39.95 Max. : 1.8291463
## open_2500.V1 med_int_2500.V1 high_int_100.V1
## Min. :-1.3794337 Min. :-2.1968789 Min. :-0.532437
## 1st Qu.:-0.5497931 1st Qu.:-0.7166467 1st Qu.:-0.532437
## Median :-0.1781892 Median : 0.2108832 Median :-0.532437
## Mean : 0.0000000 Mean : 0.0000000 Mean : 0.000000
## 3rd Qu.: 0.4072220 3rd Qu.: 0.7753690 3rd Qu.: 0.015894
## Max. : 2.7693377 Max. : 1.7719166 Max. : 4.768092
## dist_m_cafos.V1 aadt_100.V1 aadt_250.V1 x
## Min. :-2.1610758 Min. :-0.390807 Min. :-0.664113 Min. :-782614
## 1st Qu.:-0.4079546 1st Qu.:-0.390807 1st Qu.:-0.664113 1st Qu.:-763240
## Median :-0.0010458 Median :-0.390807 Median :-0.365247 Median :-756412
## Mean : 0.0000000 Mean : 0.000000 Mean : 0.000000 Mean :-756776
## 3rd Qu.: 0.3418171 3rd Qu.:-0.115098 3rd Qu.: 0.063484 3rd Qu.:-749434
## Max. : 2.2925064 Max. : 4.632471 Max. : 4.155107 Max. :-738740
## y type
## Min. :1873044 central: 1
## 1st Qu.:1890776 dist :68
## Median :1894956
## Mean :1895217
## 3rd Qu.:1899263
## Max. :1917328
cor(bc_sp_cov[,covars])
## impervious_2500 open_2500 med_int_2500 high_int_100
## impervious_2500 1.00000000 -0.8409441 0.8987424 0.48006634
## open_2500 -0.84094408 1.0000000 -0.8002841 -0.37979903
## med_int_2500 0.89874242 -0.8002841 1.0000000 0.39325577
## high_int_100 0.48006634 -0.3797990 0.3932558 1.00000000
## dist_m_cafos -0.02217726 0.3367603 -0.2416443 -0.01905821
## aadt_100 0.40911376 -0.3329512 0.2614633 0.71529092
## aadt_250 0.42742607 -0.3251319 0.2768334 0.60158329
## dist_m_cafos aadt_100 aadt_250
## impervious_2500 -0.02217726 0.4091138 0.4274261
## open_2500 0.33676028 -0.3329512 -0.3251319
## med_int_2500 -0.24164427 0.2614633 0.2768334
## high_int_100 -0.01905821 0.7152909 0.6015833
## dist_m_cafos 1.00000000 0.1015568 0.2311201
## aadt_100 0.10155682 1.0000000 0.8062597
## aadt_250 0.23112012 0.8062597 1.0000000
#' NO2, and smoke spatiotemporal predictors
#' NO2 at each site estimated using IDW
bc_st_no2 <- dplyr::select(lur_data, site_id, st_week, idw_no2) %>%
distinct() %>%
arrange(st_week) %>%
pivot_wider(id_cols = st_week,
names_from = site_id, values_from = idw_no2) %>%
as.data.frame()
rownames(bc_st_no2) <- bc_st_no2$st_week
bc_st_no2$st_week <- NULL
bc_st_no2 <- as.matrix(bc_st_no2)
no2_weeks <- rownames(bc_st_no2)
tail(no2_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
setdiff(bc_weeks, no2_weeks)
## character(0)
#' Smoke day indicator variable based on WFS paper method (see Brey and Fischer, 2016)
#' based on any smoke variable in the area using the 2 sd cut-off (area_smoke_2sd)
bc_st_smk <- dplyr::select(lur_data, site_id, st_week, area_smoke_2sd) %>%
distinct() %>%
arrange(st_week) %>%
pivot_wider(id_cols = st_week,
names_from = site_id, values_from = area_smoke_2sd) %>%
as.data.frame()
rownames(bc_st_smk) <- bc_st_smk$st_week
bc_st_smk$st_week <- NULL
bc_st_smk <- as.matrix(bc_st_smk)
smk_weeks <- rownames(bc_st_smk)
tail(smk_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
tail(bc_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
setdiff(bc_weeks, smk_weeks)
## character(0)
Create a new denver.data object.
denver.data <- createSTdata(obs = bc_obs2,
covars = bc_sp_cov,
SpatioTemporal = list(bc_st_no2 = bc_st_no2,
bc_st_smk = bc_st_smk))
D <- createDataMatrix(denver.data)
denver.data
## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 231 (observed: 231)
## No. obs: 943
##
## Only constant temporal trend,with dates:
## 2016-02-08 to 2020-11-30
##
## 13 covariate(s):
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
Next we need to combine monitoring data for NO2, BC, and PM2.5 in order to calculate the time trends.
First up: NO2
Plot the time trends for NO2 at the monitors
Original units (ppb)
no2_ts <- ggplot(data = no2_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
no2_ts
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Log-transformed NO2
no2_ts2 <- ggplot(data = no2_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
no2_ts2
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data. Use the log-transformed NO2 here
#' NO2 observations
no2_obs <- no2_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(no2 = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = no2) %>%
as.data.frame() %>%
arrange(st_week)
head(no2_obs)
head(no2_obs$st_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(no2_obs$st_week)
## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"
# NO2 SP object
no2_data2 <- read_csv(here::here("Data", "Monitor_NO2_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_no2")) %>%
filter(monitor_id %in% colnames(no2_obs))
## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
no2_coords <- do.call(rbind, st_geometry(no2_data2)) %>% as_tibble()
names(no2_coords) <- c("x", "y")
no2_sp_lonlat <- no2_data2 %>%
st_transform(crs = ll_wgs84)
no2_coords2 <- do.call(rbind, st_geometry(no2_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
no2_sp_cov <- st_set_geometry(no2_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
no2_sp_cov <- bind_cols(no2_sp_cov, no2_coords, no2_coords2) %>%
as.data.frame() %>%
distinct()
Next, BC at the central site monitor
#' BC central site concentrations
bc_cent_data <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
filter(!is.na(Arithmetic_Mean)) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(County_Code = str_sub(monitor_id, start = 3, end = 5)) %>%
filter(County_Code %in% counties) %>%
arrange(Date_Local, monitor_id) %>%
mutate(Arithmetic_Mean = ifelse(Arithmetic_Mean <= 0.005, NA, Arithmetic_Mean))
## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Pollutant_Standard = col_logical(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_logical(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
head(bc_cent_data$Date_Local)
## [1] "2016-02-08" "2016-02-09" "2016-02-10" "2016-02-11" "2016-02-12"
## [6] "2016-02-13"
tail(bc_cent_data$Date_Local)
## [1] "2020-11-25" "2020-11-26" "2020-11-27" "2020-11-28" "2020-11-29"
## [6] "2020-11-30"
#' Add an identifier to differentiate these measurements from collocated pollutants at the same site
bc_cent_data$monitor_id <- paste0(bc_cent_data$monitor_id, "_bc")
unique(bc_cent_data$monitor_id)
## [1] "080310027_bc"
Plot the time trends for BC at the central site monitor
Original units (ug/m3)
bc_ts <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
#facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
bc_ts
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Log-transformed BC
bc_ts2 <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
#facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
bc_ts2
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data
Use log-transformed data here
#' BC observations
bc_cent_obs <- bc_cent_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(bc = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = bc) %>%
as.data.frame() %>%
arrange(st_week)
head(bc_cent_obs)
head(bc_cent_obs$st_week)
## [1] "2016-02-08" "2016-02-15" "2016-02-22" "2016-02-29" "2016-03-07"
## [6] "2016-03-14"
tail(bc_cent_obs$st_week)
## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"
# BC Central Site SP object
bc_cent_data2 <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_bc")) %>%
filter(monitor_id %in% colnames(bc_cent_obs))
## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Pollutant_Standard = col_logical(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_logical(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
bc_cent_coords <- do.call(rbind, st_geometry(bc_cent_data2)) %>% as_tibble()
names(bc_cent_coords) <- c("x", "y")
bc_cent_sp_lonlat <- bc_cent_data2 %>%
st_transform(crs = ll_wgs84)
bc_cent_coords2 <- do.call(rbind, st_geometry(bc_cent_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
bc_cent_sp_cov <- st_set_geometry(bc_cent_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
bc_cent_sp_cov <- bind_cols(bc_cent_sp_cov, bc_cent_coords, bc_cent_coords2) %>%
as.data.frame() %>%
distinct()
Lastly, PM2.5 from select monitors with long-term records
Plot the time trends for PM2.5 at the monitors
Original units (ug/m3)
pm_ts <- ggplot(data = pm_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
pm_ts
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Log-transformed PM2.5
pm_ts2 <- ggplot(data = pm_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
pm_ts2
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data
Use the log-transformed data here
#' PM2.5 observations
pm_obs <- pm_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(pm = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = pm) %>%
as.data.frame() %>%
arrange(st_week)
head(pm_obs)
head(pm_obs$st_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(pm_obs$st_week)
## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"
# PM SP object
pm_data2 <- read_csv(here::here("Data", "Monitor_PM_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_pm")) %>%
filter(monitor_id %in% colnames(pm_obs))
## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_double(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
pm_coords <- do.call(rbind, st_geometry(pm_data2)) %>% as_tibble()
names(pm_coords) <- c("x", "y")
pm_sp_lonlat <- pm_data2 %>%
st_transform(crs = ll_wgs84)
pm_coords2 <- do.call(rbind, st_geometry(pm_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
pm_sp_cov <- st_set_geometry(pm_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
pm_sp_cov <- bind_cols(pm_sp_cov, pm_coords, pm_coords2) %>%
as.data.frame() %>%
distinct()
Combine the observations and the spatial covariates Scale the pollutant measurements
#' Join log-transformed observations
pol_obs <- data.frame(st_week = rownames(bc_obs2)) %>%
left_join(no2_obs, by = "st_week") %>%
left_join(bc_cent_obs, by = "st_week") %>%
left_join(pm_obs, by = "st_week")
glimpse(pol_obs)
## Rows: 627
## Columns: 18
## $ st_week <chr> "2008-12-29", "2009-01-05", "2009-01-12", "2009-01-19…
## $ `080013001_no2` <dbl> 2.4313422, 2.3337153, 3.0953820, 2.4960668, 2.6445846…
## $ `080310002_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310026_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_bc` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080010006_pm` <dbl> 1.394118, 1.247789, 1.934844, 2.298456, 1.598321, 1.7…
## $ `080010008_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080050005_pm` <dbl> 1.2490759, 0.6817687, 1.3697744, 1.7730816, 1.1939225…
## $ `080130003_pm` <dbl> 0.9679299, 0.8243293, 2.1047287, 1.7212142, 1.9059908…
## $ `080130012_pm` <dbl> 1.0919008, 0.8908546, 1.2511276, 1.2703657, 0.8047190…
## $ `080310002_pm` <dbl> 1.335209, 1.188381, 1.767997, 2.027844, 1.486213, 1.7…
## $ `080310023_pm` <dbl> 1.661197, 1.513911, 1.919262, 1.996528, 1.909773, 1.7…
## $ `080310025_pm` <dbl> 1.1537863, 0.9359011, 1.3262687, 1.9969858, 1.6521301…
## $ `080310026_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
rownames(pol_obs) <- pol_obs$st_week
pol_obs$st_week <- NULL
pol_obs <- as.matrix(pol_obs)
#' Scale the measurements to avoid issues where units are different
#' Remember! These have been log-transformed before scaling
#' Scale the data by columns (monitors)
rnames <- rownames(pol_obs)
cnames <- colnames(pol_obs)
pol_obs <- apply(pol_obs, 2, scale)
rownames(pol_obs) <- rnames
colnames(pol_obs) <- cnames
head(rownames(pol_obs))
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(rownames(pol_obs))
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
head(colnames(pol_obs))
## [1] "080013001_no2" "080310002_no2" "080310026_no2" "080310027_no2"
## [5] "080310028_no2" "080310027_bc"
class(pol_obs)
## [1] "matrix" "array"
dim(pol_obs)
## [1] 627 17
#' SP variables
pol_sp_cov <- bind_rows(no2_sp_cov, bc_cent_sp_cov, pm_sp_cov)
pol_sp_cov
Create the data object
pol_STdata <- createSTdata(obs = pol_obs,
covars = pol_sp_cov)
print(pol_STdata)
## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 624 (observed: 624)
## No. obs: 6574
##
## Only constant temporal trend,with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.
The cross validation results for generating the time trends suggest 2 basis functions is the way to go, but let’s see if that holds up with our data
D_pol <- createDataMatrix(pol_STdata)
#D_pol
n_years <- length(2009:2020)
n_years*4
## [1] 48
n_years*8
## [1] 96
#' Trying 4 df per year
pol.SVD.cv.4py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 4)
print(pol.SVD.cv.4py)
## Result of SVDsmoothCV, average of CV-statistics:
## MSE R2 AIC BIC
## n.basis.0 0.9975662 0.0000000 0.9987807 5.097602
## n.basis.1 0.8016075 0.1964669 -96.6267210 -88.429078
## n.basis.2 0.6359015 0.3625417 -225.5766119 -213.280148
## n.basis.3 0.6299681 0.3684906 -229.4640525 -213.068767
## n.basis.4 0.6089310 0.3896015 -240.3947203 -219.900613
plot(pol.SVD.cv.4py)
#' Trying 8 df per year
pol.SVD.cv.8py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 8)
print(pol.SVD.cv.8py)
## Result of SVDsmoothCV, average of CV-statistics:
## MSE R2 AIC BIC
## n.basis.0 0.9975662 0.0000000 0.9987807 5.097602
## n.basis.1 0.7615784 0.2365880 -121.3718836 -113.174241
## n.basis.2 0.5808908 0.4176809 -266.8631552 -254.566691
## n.basis.3 0.5724708 0.4261239 -272.6314965 -256.236211
## n.basis.4 0.5602312 0.4384013 -279.8920086 -259.397901
plot(pol.SVD.cv.8py)
I tried using 1 and 2 basis functions for the time trends. I also checked out 4 and 8 df per year
First, see what the time trends for the monitoring data look like with 1 basis function.
4 df per year
pol_STdata1 <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 4)
## Replacing existing trend.
pol_STdata1
## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 624 (observed: 624)
## No. obs: 6574
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.
head(pol_STdata1$trend)
tail(pol_STdata1$trend)
plot(pol_STdata1$trend$date, pol_STdata1$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1$trend$date, pol_STdata1$atrend$V1, col = 1)
8 df per year
pol_STdata1a <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 8)
## Replacing existing trend.
pol_STdata1a
## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 624 (observed: 624)
## No. obs: 6574
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.
head(pol_STdata1a$trend)
tail(pol_STdata1a$trend)
plot(pol_STdata1a$trend$date, pol_STdata1a$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1a$trend$date, pol_STdata1a$atrend$V1, col = 1)
Next, let’s see what 2 basis functions look like:
4 df per year
pol_STdata2 <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 4)
## Replacing existing trend.
pol_STdata2
## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 624 (observed: 624)
## No. obs: 6574
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.
head(pol_STdata2$trend)
tail(pol_STdata2$trend)
plot(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2)
8 df per year
pol_STdata2a <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 8)
## Replacing existing trend.
pol_STdata2a
## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 624 (observed: 624)
## No. obs: 6574
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.
head(pol_STdata2a$trend)
tail(pol_STdata2a$trend)
plot(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2)
The next step is to update the trend data for the denver.data object and see if these trend data fit our observations.
Comparing the plots using one basis function vs. two basis functions, it looks like using two basis functions might be better than 1. There was some residual autocorrelation in the central site monitor data when using just one basis function (see the fourth plot in the first series of plots below). Using 8 degrees of freedom did not help.
Using two basis functions helped somewhat to reduce this autocorrelation, but it still persists. There are still some noticable trends in the residuals for the central monitoring site
First update the trend data, then plot the trends for three distributed sites and the central BC monitoring site.
denver.data2.1 <- denver.data
denver.data2.1$trend <- pol_STdata1$trend
print(denver.data2.1)
## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 624 (observed: 231)
## No. obs: 943
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 13 covariate(s):
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.1, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_1")
plot(denver.data2.1, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_10")
plot(denver.data2.1, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_43")
plot(denver.data2.1, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.1, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_53")
plot(denver.data2.1, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.1, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="central")
plot(denver.data2.1, "pacf", ID="central")
denver.data2.1a <- denver.data
denver.data2.1a$trend <- pol_STdata1a$trend
print(denver.data2.1a)
## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 624 (observed: 231)
## No. obs: 943
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 13 covariate(s):
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.1a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_1")
plot(denver.data2.1a, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_10")
plot(denver.data2.1a, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_43")
plot(denver.data2.1a, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.1a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_53")
plot(denver.data2.1a, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.1a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="central")
plot(denver.data2.1a, "pacf", ID="central")
denver.data2.2 <- denver.data
denver.data2.2$trend <- pol_STdata2$trend
print(denver.data2.2)
## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 624 (observed: 231)
## No. obs: 943
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 13 covariate(s):
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.2, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_1")
plot(denver.data2.2, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_10")
plot(denver.data2.2, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_10")
plot(denver.data2.2, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.2, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_43")
plot(denver.data2.2, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.2, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_53")
plot(denver.data2.2, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.2, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="central")
plot(denver.data2.2, "pacf", ID="central")
denver.data2.2a <- denver.data
denver.data2.2a$trend <- pol_STdata2a$trend
print(denver.data2.2a)
## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 624 (observed: 231)
## No. obs: 943
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 13 covariate(s):
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.2a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_1")
plot(denver.data2.2a, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_10")
plot(denver.data2.2a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_10")
plot(denver.data2.2a, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.2a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_43")
plot(denver.data2.2a, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.2a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_53")
plot(denver.data2.2a, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.2a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="central")
plot(denver.data2.2a, "pacf", ID="central")
Building on the results of the previous modeling efforts, I’m going to update “Model 4.2”, which used an ’exp covariance structure for the error term and iid for the beta fields.
For this version of the model, use iid for cov.beta (beta, beta1) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.
denver.data.A <- denver.data2.1
names(denver.data.A$covars)
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
LUR.A <- list(covar_fun, covar_fun)
cov.beta.A <- list(covf = c("iid", "iid"), nugget = T)
cov.nu.A <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.A <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")
denver.model.A <- createSTmodel(denver.data.A, LUR = LUR.A,
ST = c("bc_st_no2", "bc_st_smk"),
cov.beta = cov.beta.A, cov.nu = cov.nu.A,
locations = locations.A)
## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.
denver.model.A
## STmodel-object with:
## No. locations: 69 (observed: 64)
## No. time points: 624 (observed: 231)
## No. obs: 943
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## Models for the beta-fields are:
## $const
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 +
## dist_m_cafos + aadt_100 + aadt_250
##
## $V1
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 +
## dist_m_cafos + aadt_100 + aadt_250
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## Covariance model for the beta-field(s):
## Covariance type(s): iid, iid
## Nugget: Yes, Yes
## Covariance model for the nu-field(s):
## Covariance type: exp
## Nugget: ~1
## Random effect: No
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
set.seed(1000)
names <- loglikeSTnames(denver.model.A, all=FALSE)
names
## [1] "log.nugget.const.iid" "log.nugget.V1.iid"
## [3] "nu.log.range.exp" "nu.log.sill.exp"
## [5] "nu.log.nugget.(Intercept).exp"
x.init.A <- cbind(c(0, 0, 0, 0, 0),
c(-1, -1, 4, -5, -1),
c(-5, -5, 4, -1, -5),
c(-1, -1, 8, -5, -1),
c(-5, -5, 8, -1, -5),
c(-1, -1, 12, -5, -1),
c(-5, -5, 12, -1, -5),
c(-7, -7, 12, -3, -3))
rownames(x.init.A) <- loglikeSTnames(denver.model.A, all=FALSE)
x.init.A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid 0 -1 -5 -1 -5 -1 -5 -7
## log.nugget.V1.iid 0 -1 -5 -1 -5 -1 -5 -7
## nu.log.range.exp 0 4 4 8 8 12 12 12
## nu.log.sill.exp 0 -5 -1 -5 -1 -5 -1 -3
## nu.log.nugget.(Intercept).exp 0 -1 -5 -1 -5 -1 -5 -3
#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.A <- estimate.STmodel(denver.model.A, x.init.A)
## Optimisation using starting value 1/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 428.12 |proj g|= 15
## At iterate 10 f = -1187.6 |proj g|= 1.5251
##
## iterations 18
## function evaluations 36
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0128187
## final function value -1187.79
##
## F = -1187.79
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1187.786129
## converged
## Optimisation using starting value 2/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -299.6 |proj g|= 14
## At iterate 10 f = -1184.3 |proj g|= 0.423
## ys=-1.276e-02 -gs= 2.763e-02, BFGS update SKIPPED
## At iterate 20 f = -1184.9 |proj g|= 11.462
## At iterate 30 f = -1187.8 |proj g|= 0.19815
## ys=-8.732e-01 -gs= 1.086e+00, BFGS update SKIPPED
## At iterate 40 f = -1403.6 |proj g|= 0.18215
##
## iterations 43
## function evaluations 76
## segments explored during Cauchy searches 46
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00508496
## final function value -1403.62
##
## F = -1403.62
## final value -1403.623292
## converged
## Optimisation using starting value 3/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -421.29 |proj g|= 14
## At iterate 10 f = -1183.3 |proj g|= 1.6336
## At iterate 20 f = -1187.7 |proj g|= 0.6575
## ys=-7.955e+01 -gs= 7.565e-01, BFGS update SKIPPED
## At iterate 30 f = -1401 |proj g|= 10.6
## At iterate 40 f = -1409.4 |proj g|= 0.36759
##
## iterations 48
## function evaluations 87
## segments explored during Cauchy searches 49
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0252493
## final function value -1409.4
##
## F = -1409.4
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1409.403118
## converged
## Optimisation using starting value 4/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -299.84 |proj g|= 14
## At iterate 10 f = -1400.2 |proj g|= 7.974
## At iterate 20 f = -1405 |proj g|= 7.8782
## At iterate 30 f = -1407.6 |proj g|= 1.8038
## At iterate 40 f = -1409.3 |proj g|= 2.6564
## At iterate 50 f = -1409.4 |proj g|= 0.00075243
##
## iterations 50
## function evaluations 67
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00075243
## final function value -1409.4
##
## F = -1409.4
## final value -1409.403111
## converged
## Optimisation using starting value 5/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -495.77 |proj g|= 14
## At iterate 10 f = -1409.3 |proj g|= 0.80461
## At iterate 20 f = -1409.4 |proj g|= 0.15163
##
## iterations 24
## function evaluations 36
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0995913
## final function value -1409.4
##
## F = -1409.4
## final value -1409.403060
## converged
## Optimisation using starting value 6/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -302.63 |proj g|= 14
## At iterate 10 f = -1408.9 |proj g|= 8.1359
## At iterate 20 f = -1409.4 |proj g|= 0.0029418
##
## iterations 21
## function evaluations 29
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.010386
## final function value -1409.4
##
## F = -1409.4
## final value -1409.403107
## converged
## Optimisation using starting value 7/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1152.5 |proj g|= 14
## At iterate 10 f = -1409.4 |proj g|= 0.33301
## At iterate 20 f = -1409.4 |proj g|= 0.025957
##
## iterations 20
## function evaluations 40
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.025957
## final function value -1409.4
##
## F = -1409.4
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1409.403103
## converged
## Optimisation using starting value 8/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1118.9 |proj g|= 12
## At iterate 10 f = -1409.4 |proj g|= 0.13337
##
## iterations 15
## function evaluations 23
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.127558
## final function value -1409.4
##
## F = -1409.4
## final value -1409.402987
## converged
print(est.denver.model.A)
## Optimisation for STmodel with 8 starting points.
## Results: 6 converged, 2 not converged, 0 failed.
## Best result for starting point 3, optimisation has converged
##
## No fixed parameters.
##
## Estimated parameters for all starting point(s):
## [,1] [,2] [,3]
## gamma.bc_st_no2 0.0190177471 0.0139611091 0.013524343
## gamma.bc_st_smk -0.0426958408 -0.0160214985 -0.017276051
## alpha.const.(Intercept) -0.1503355265 0.0025614497 0.008365957
## alpha.const.impervious_2500 0.0376285789 0.0297836107 0.030578606
## alpha.const.open_2500 -0.0076291505 -0.0104722683 -0.010608880
## alpha.const.med_int_2500 -0.0164412434 -0.0205190931 -0.020872498
## alpha.const.high_int_100 0.0113010906 0.0123640421 0.010811603
## alpha.const.dist_m_cafos -0.0373287719 -0.0354396137 -0.035629090
## alpha.const.aadt_100 0.0322796488 0.0292214710 0.029130554
## alpha.const.aadt_250 0.0164520345 0.0148726789 0.014634801
## alpha.V1.(Intercept) 0.1874787211 0.2202093765 0.225017101
## alpha.V1.impervious_2500 0.0117787671 -0.0089070999 -0.009956409
## alpha.V1.open_2500 0.0118058417 0.0114616390 0.011181075
## alpha.V1.med_int_2500 0.0090074275 0.0130834301 0.011830049
## alpha.V1.high_int_100 0.0180754356 0.0178398056 0.018430627
## alpha.V1.dist_m_cafos -0.0002816461 -0.0004239591 -0.002006719
## alpha.V1.aadt_100 -0.0323430204 -0.0319474609 -0.030466145
## alpha.V1.aadt_250 -0.0150887087 -0.0127119439 -0.010031267
## log.nugget.const.iid -14.7001121803 -14.9989138692 -7.558680048
## log.nugget.V1.iid -15.0000000000 -14.9994643477 -7.132387053
## nu.log.range.exp 0.0000000000 12.5059102289 12.520026954
## nu.log.sill.exp -4.3239982804 -3.0584488703 -3.032854210
## nu.log.nugget.(Intercept).exp -4.1033971401 -4.6498599256 -4.759856021
## [,4] [,5] [,6]
## gamma.bc_st_no2 0.013525741 0.013529300 0.013525486
## gamma.bc_st_smk -0.017271850 -0.017263869 -0.017272711
## alpha.const.(Intercept) 0.008340578 0.008261481 0.008346703
## alpha.const.impervious_2500 0.030578311 0.030574416 0.030578433
## alpha.const.open_2500 -0.010609062 -0.010610972 -0.010609050
## alpha.const.med_int_2500 -0.020872571 -0.020871243 -0.020872644
## alpha.const.high_int_100 0.010812892 0.010811989 0.010812616
## alpha.const.dist_m_cafos -0.035630328 -0.035631806 -0.035630314
## alpha.const.aadt_100 0.029131348 0.029132930 0.029131662
## alpha.const.aadt_250 0.014630979 0.014625332 0.014630607
## alpha.V1.(Intercept) 0.225008013 0.224984562 0.225008189
## alpha.V1.impervious_2500 -0.009955584 -0.009961123 -0.009956377
## alpha.V1.open_2500 0.011183188 0.011185787 0.011183151
## alpha.V1.med_int_2500 0.011827663 0.011826552 0.011827449
## alpha.V1.high_int_100 0.018433956 0.018435501 0.018434386
## alpha.V1.dist_m_cafos -0.002006304 -0.002010938 -0.002006642
## alpha.V1.aadt_100 -0.030462478 -0.030452721 -0.030461775
## alpha.V1.aadt_250 -0.010032886 -0.010031414 -0.010032386
## log.nugget.const.iid -7.560464968 -7.560107182 -7.560275086
## log.nugget.V1.iid -7.131144030 -7.127587858 -7.130684472
## nu.log.range.exp 12.519691760 12.517936372 12.519742836
## nu.log.sill.exp -3.033134958 -3.033381695 -3.033081239
## nu.log.nugget.(Intercept).exp -4.759862263 -4.759736916 -4.759927367
## [,7] [,8]
## gamma.bc_st_no2 0.013523845 0.013523998
## gamma.bc_st_smk -0.017273875 -0.017269878
## alpha.const.(Intercept) 0.008380116 0.008390079
## alpha.const.impervious_2500 0.030583778 0.030585150
## alpha.const.open_2500 -0.010607883 -0.010605271
## alpha.const.med_int_2500 -0.020873883 -0.020875042
## alpha.const.high_int_100 0.010808899 0.010821491
## alpha.const.dist_m_cafos -0.035630538 -0.035630373
## alpha.const.aadt_100 0.029132549 0.029126451
## alpha.const.aadt_250 0.014630832 0.014637730
## alpha.V1.(Intercept) 0.225021583 0.225035682
## alpha.V1.impervious_2500 -0.009955880 -0.009932152
## alpha.V1.open_2500 0.011181831 0.011183267
## alpha.V1.med_int_2500 0.011824019 0.011826210
## alpha.V1.high_int_100 0.018438465 0.018434686
## alpha.V1.dist_m_cafos -0.002005353 -0.001990822
## alpha.V1.aadt_100 -0.030465340 -0.030482658
## alpha.V1.aadt_250 -0.010026484 -0.010046108
## log.nugget.const.iid -7.557186151 -7.567336821
## log.nugget.V1.iid -7.127421757 -7.140153420
## nu.log.range.exp 12.520867210 12.522302555
## nu.log.sill.exp -3.032945560 -3.033012579
## nu.log.nugget.(Intercept).exp -4.759776099 -4.758603669
##
## Function value(s):
## [1] 1187.786 1403.623 1409.403 1409.403 1409.403 1409.403 1409.403 1409.403
Define the CV groups
set.seed(1000)
site_idsA <- unique(denver.model.A$obs$ID[which(str_detect(denver.model.A$obs$ID, "d_") == T)])
site_idsA
## [1] "d_1" "d_18" "d_22" "d_25" "d_28" "d_31" "d_33" "d_36" "d_39" "d_4"
## [11] "d_7" "d_11" "d_12" "d_19" "d_23" "d_26" "d_29" "d_32" "d_35" "d_38"
## [21] "d_40" "d_8" "d_42" "d_44" "d_46" "d_50" "d_52" "d_54" "d_56" "d_58"
## [31] "d_61" "d_63" "d_65" "d_68" "d_15" "d_41" "d_43" "d_45" "d_47" "d_51"
## [41] "d_53" "d_67" "d_49" "d_14" "d_48" "d_37" "d_55" "d_57" "d_64" "d_66"
## [51] "d_16" "d_2" "d_5" "d_60" "d_62" "d_34" "d_10" "d_13" "d_21" "d_3"
## [61] "d_6" "d_17" "d_20"
Ind.cv.A <- createCV(denver.model.A, groups = 10, min.dist = .1,
subset = site_idsA)
ID.cv.A <- sapply(split(denver.model.A$obs$ID, Ind.cv.A), unique)
print(sapply(ID.cv.A, length))
## 0 1 2 3 4 5 6 7 8 9 10
## 1 7 7 7 6 6 6 6 6 6 6
table(Ind.cv.A)
## Ind.cv.A
## 0 1 2 3 4 5 6 7 8 9 10
## 231 66 76 82 67 57 85 72 78 54 75
I.col.A <- apply(sapply(ID.cv.A, function(x) denver.model.A$locations$ID%in% x), 1,
function(x) if(sum(x)==1) which(x) else 0)
names(I.col.A) <- denver.model.A$locations$ID
print(I.col.A)
## central d_1 d_10 d_11 d_12 d_13 d_14 d_15 d_16 d_17
## 1 6 2 3 4 4 4 9 4 10
## d_18 d_19 d_2 d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 8 5 11 3 2 11 4 9 11 9
## d_29 d_3 d_31 d_32 d_33 d_34 d_35 d_36 d_37 d_38
## 2 4 6 8 5 5 10 8 3 8
## d_39 d_4 d_40 d_41 d_42 d_43 d_44 d_45 d_46 d_47
## 6 7 7 3 7 8 11 3 9 5
## d_48 d_49 d_5 d_50 d_51 d_52 d_53 d_54 d_55 d_56
## 9 4 6 9 10 3 7 5 7 6
## d_57 d_58 d_6 d_60 d_61 d_62 d_63 d_64 d_65 d_66
## 2 2 11 5 2 11 10 3 8 6
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 7 10 2 10 0 0 0 0 0
par(mfrow=c(1,1))
plot(denver.model.A$locations$long,
denver.model.A$locations$lat,
pch=23+floor(I.col.A/max(I.col.A)+.5), bg=I.col.A,
xlab="Longitude", ylab="Latitude")
map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)
## [[1]]
## NULL
ID starting conditions using previous model without CV:
#' ID starting conditions using previous model without CV:
x.init.A.cv <- coef(est.denver.model.A, pars="cov")[,c("par","init")]
x.init.A.cv
Run the model with cross validation
est.denver.A.cv <- estimateCV(denver.model.A, x.init.A.cv, Ind.cv.A)
##
## ***************************
## Estimation of CV-set 1/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1279.4 |proj g|= 13.463
## At iterate 10 f = -1279.9 |proj g|= 0.17604
##
## iterations 16
## function evaluations 24
## segments explored during Cauchy searches 16
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00869475
## final function value -1279.92
##
## F = -1279.92
## final value -1279.916549
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -391.08 |proj g|= 14
## At iterate 10 f = -1086 |proj g|= 2.3228
## ys=-4.014e-01 -gs= 1.977e-01, BFGS update SKIPPED
## At iterate 20 f = -1258.2 |proj g|= 19.721
## At iterate 30 f = -1275.4 |proj g|= 0.0050604
##
## iterations 30
## function evaluations 49
## segments explored during Cauchy searches 31
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00506041
## final function value -1275.42
##
## F = -1275.42
## final value -1275.424691
## converged
##
##
## ***************************
## Estimation of CV-set 2/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1280.8 |proj g|= 1.2959
## At iterate 10 f = -1280.9 |proj g|= 0.40948
##
## iterations 15
## function evaluations 31
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0212405
## final function value -1280.93
##
## F = -1280.93
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1280.932738
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -387.84 |proj g|= 14
## At iterate 10 f = -1086.5 |proj g|= 2.4518
## At iterate 20 f = -1090.4 |proj g|= 3.3494
## At iterate 30 f = -1272.2 |proj g|= 3.3665
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1275.520364
## stopped after 37 iterations
##
##
## ***************************
## Estimation of CV-set 3/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1285.1 |proj g|= 10.24
##
## iterations 9
## function evaluations 18
## segments explored during Cauchy searches 11
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0629467
## final function value -1285.94
##
## F = -1285.94
## final value -1285.938947
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -378.99 |proj g|= 14
## At iterate 10 f = -1080.4 |proj g|= 13.47
## ys=-3.594e-01 -gs= 6.181e-01, BFGS update SKIPPED
## ys=-2.044e+00 -gs= 2.394e+00, BFGS update SKIPPED
## At iterate 20 f = -1196 |proj g|= 18.249
## At iterate 30 f = -1276.8 |proj g|= 12.846
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1278.544393
## stopped after 37 iterations
##
##
## ***************************
## Estimation of CV-set 4/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1273.2 |proj g|= 16.881
## At iterate 10 f = -1274 |proj g|= 0.15717
##
## iterations 17
## function evaluations 34
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0382929
## final function value -1274.04
##
## F = -1274.04
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1274.042573
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -390.46 |proj g|= 14
## At iterate 10 f = -1083 |proj g|= 2.3331
## At iterate 20 f = -1087.3 |proj g|= 0.31284
## ys=-6.049e-01 -gs= 1.042e+00, BFGS update SKIPPED
## At iterate 30 f = -1147.2 |proj g|= 20.9
## At iterate 40 f = -1269 |proj g|= 10.561
##
## iterations 45
## function evaluations 65
## segments explored during Cauchy searches 46
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00406416
## final function value -1270.08
##
## F = -1270.08
## final value -1270.080154
## converged
##
##
## ***************************
## Estimation of CV-set 5/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1300.8 |proj g|= 9.1053
## At iterate 10 f = -1301.2 |proj g|= 0.080753
##
## iterations 17
## function evaluations 40
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0469096
## final function value -1301.24
##
## F = -1301.24
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1301.239089
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -395.42 |proj g|= 14
## At iterate 10 f = -1100.6 |proj g|= 1.8267
## At iterate 20 f = -1105 |proj g|= 0.36357
## ys=-2.770e+00 -gs= 4.355e-02, BFGS update SKIPPED
## At iterate 30 f = -1292.1 |proj g|= 19.885
##
## iterations 38
## function evaluations 50
## segments explored during Cauchy searches 39
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0040607
## final function value -1297.73
##
## F = -1297.73
## final value -1297.728676
## converged
##
##
## ***************************
## Estimation of CV-set 6/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1273.3 |proj g|= 8.8806
## At iterate 10 f = -1273.6 |proj g|= 0.035777
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 16
## function evaluations 42
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00197049
## final function value -1273.62
##
## F = -1273.62
## final value -1273.617044
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -382.95 |proj g|= 14
## At iterate 10 f = -1065 |proj g|= 3.3156
## At iterate 20 f = -1069.2 |proj g|= 0.27056
## At iterate 30 f = -1267.7 |proj g|= 5.0383
##
## iterations 36
## function evaluations 61
## segments explored during Cauchy searches 37
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0236423
## final function value -1268.61
##
## F = -1268.61
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1268.613291
## converged
##
##
## ***************************
## Estimation of CV-set 7/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1281.2 |proj g|= 7.644
## At iterate 10 f = -1282.8 |proj g|= 0.13879
## At iterate 20 f = -1282.8 |proj g|= 0.039849
##
## iterations 20
## function evaluations 37
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0398489
## final function value -1282.79
##
## F = -1282.79
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1282.788383
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -388.62 |proj g|= 14
## At iterate 10 f = -1080.1 |proj g|= 2.3357
## ys=-4.483e-02 -gs= 2.380e-01, BFGS update SKIPPED
## ys=-1.945e+02 -gs= 6.398e+00, BFGS update SKIPPED
## At iterate 20 f = -1260.6 |proj g|= 8.8126
## At iterate 30 f = -1280 |proj g|= 0.19367
##
## iterations 32
## function evaluations 43
## segments explored during Cauchy searches 33
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.116086
## final function value -1280.01
##
## F = -1280.01
## final value -1280.007071
## converged
##
##
## ***************************
## Estimation of CV-set 8/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1263.8 |proj g|= 9.6277
## At iterate 10 f = -1265 |proj g|= 0.18815
## At iterate 20 f = -1265 |proj g|= 0.015515
##
## iterations 22
## function evaluations 37
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0126596
## final function value -1265.02
##
## F = -1265.02
## final value -1265.019277
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -385.81 |proj g|= 14
## At iterate 10 f = -1070.3 |proj g|= 1.7684
## At iterate 20 f = -1074.7 |proj g|= 0.64859
## ys=-1.518e-01 -gs= 1.371e-01, BFGS update SKIPPED
## ys=-7.603e+01 -gs= 3.615e+00, BFGS update SKIPPED
## At iterate 30 f = -1177.2 |proj g|= 19.527
## At iterate 40 f = -1262 |proj g|= 0.004444
##
## iterations 41
## function evaluations 64
## segments explored during Cauchy searches 42
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0044462
## final function value -1261.98
##
## F = -1261.98
## final value -1261.975372
## converged
##
##
## ***************************
## Estimation of CV-set 9/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1305.7 |proj g|= 8.8586
## At iterate 10 f = -1306 |proj g|= 0.26059
##
## iterations 13
## function evaluations 26
## segments explored during Cauchy searches 13
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.051015
## final function value -1306.04
##
## F = -1306.04
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1306.044405
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -396.77 |proj g|= 14
## At iterate 10 f = -1105.6 |proj g|= 2.2523
## At iterate 20 f = -1109.9 |proj g|= 0.25648
## ys=-6.297e-01 -gs= 1.942e+00, BFGS update SKIPPED
## At iterate 30 f = -1302.6 |proj g|= 0.55208
##
## iterations 35
## function evaluations 52
## segments explored during Cauchy searches 36
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00787603
## final function value -1302.6
##
## F = -1302.6
## final value -1302.604310
## converged
##
##
## ***************************
## Estimation of CV-set 10/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1292.3 |proj g|= 10.24
## At iterate 10 f = -1293.8 |proj g|= 0.61724
## At iterate 20 f = -1293.8 |proj g|= 0.020958
##
## iterations 20
## function evaluations 41
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0209579
## final function value -1293.76
##
## F = -1293.76
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1293.755426
## converged
##
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -387.69 |proj g|= 14
## At iterate 10 f = -1081.4 |proj g|= 2.9417
## At iterate 20 f = -1085.5 |proj g|= 0.22294
##
## iterations 23
## function evaluations 37
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0124233
## final function value -1085.47
##
## F = -1085.47
## final value -1085.473201
## converged
##
print(est.denver.A.cv)
## Cross-validation parameter estimation for STmodel
## with 10 CV-groups and 2 starting points.
## Results: 10 converged, 0 not converged.
##
## No fixed parameters.
##
## Estimated function values and convergence info:
## value convergence conv eigen.min eigen.all.min
## 1 1279.917 TRUE TRUE 1.8351328 NA
## 2 1280.933 TRUE TRUE 2.7925040 NA
## 3 1285.939 TRUE TRUE 2.9117065 NA
## 4 1274.043 TRUE TRUE 1.4687175 NA
## 5 1301.239 TRUE TRUE 1.3544770 NA
## 6 1273.617 TRUE TRUE 1.7810323 NA
## 7 1282.788 TRUE TRUE 1.3687804 NA
## 8 1265.019 TRUE TRUE 0.2122064 NA
## 9 1306.044 TRUE TRUE 1.8032551 NA
## 10 1293.755 TRUE TRUE 3.1516006 NA
par(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.A, pars="all"),
plotCI((1:length(par))+.3, par, uiw=1.96*sd,
col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.A.cv, "all", boxwex=.4, col="grey", add=TRUE)
#' Save the results as an .rdata object
save(denver.data.A, denver.model.A, est.denver.model.A, est.denver.A.cv,
file = here::here("Results", "Denver_ST_Model_A.rdata"))
Making predictions using the CV model. Printing out the CV summary statistics as well
pred.A.cv <- predictCV(denver.model.A, est.denver.A.cv, LTA = T)
pred.A.cv.log <- predictCV(denver.model.A, est.denver.A.cv, LTA = T,
transform="unbiased")
head(pred.A.cv$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13
## 2008-12-29 NA -0.42509020 -0.29919410 -0.253778802 -0.29997252 -0.28480778
## 2009-01-05 NA -0.33745401 -0.22355879 -0.176026184 -0.22340549 -0.20997248
## 2009-01-12 NA -0.16860340 -0.06722598 -0.017699541 -0.07023096 -0.05849849
## 2009-01-19 NA -0.19323774 -0.10221236 -0.050581637 -0.09906798 -0.08896623
## 2009-01-26 NA -0.12695671 -0.04638523 0.007117754 -0.04189279 -0.03331327
## 2009-02-02 NA -0.03737667 0.03344767 0.088619327 0.03771543 0.04491357
## d_14 d_15 d_16 d_17 d_18
## 2008-12-29 -0.46606666 -0.388658381 0.2034593 -0.398988727 -0.4093148
## 2009-01-05 -0.38294825 -0.298634900 0.2514619 -0.309999522 -0.3254962
## 2009-01-12 -0.22334034 -0.124053289 0.3765865 -0.141111593 -0.1605016
## 2009-01-19 -0.24600816 -0.150822997 0.3208514 -0.162800490 -0.1885920
## 2009-01-26 -0.18307413 -0.082883369 0.3529177 -0.095162411 -0.1256528
## 2009-02-02 -0.09823997 0.009405756 0.4097405 -0.004809943 -0.0391435
## d_19 d_2 d_20 d_21 d_22
## 2008-12-29 -0.26739332 0.2419826 -0.233183779 -0.38157877 -0.43694213
## 2009-01-05 -0.18619766 0.2934278 -0.159571444 -0.29682928 -0.34550831
## 2009-01-12 -0.01872363 0.4280098 -0.005310516 -0.13154643 -0.17165784
## 2009-01-19 -0.05582689 0.3676511 -0.042091373 -0.15795028 -0.19436063
## 2009-01-26 0.00415006 0.4019541 0.011968596 -0.09411154 -0.12490642
## 2009-02-02 0.08976834 0.4629670 0.090167521 -0.00700837 -0.03199512
## d_23 d_25 d_26 d_28 d_29
## 2008-12-29 -0.03727008 -0.322228814 -0.31186712 -0.40101437 -0.5158686
## 2009-01-05 0.02486235 -0.246537633 -0.23289928 -0.32203488 -0.4293773
## 2009-01-12 0.16386225 -0.086030188 -0.07129028 -0.15829835 -0.2623840
## 2009-01-19 0.12143267 -0.126296122 -0.10573185 -0.19546780 -0.2871477
## 2009-01-26 0.16591946 -0.070954973 -0.04723558 -0.13723614 -0.2217779
## 2009-02-02 0.23401340 0.009901469 0.03573177 -0.05375665 -0.1332853
## d_3 d_31 d_32 d_33 d_34
## 2008-12-29 -0.30794846 -0.13373488 -0.5138657 -0.43435799 -0.236532203
## 2009-01-05 -0.22897077 -0.06982084 -0.4273476 -0.34848710 -0.162280634
## 2009-01-12 -0.07342900 0.07573486 -0.2597022 -0.17642206 -0.001625635
## 2009-01-19 -0.09999600 0.02876220 -0.2852506 -0.20912282 -0.045267901
## 2009-01-26 -0.04070177 0.07419080 -0.2199385 -0.14503622 0.008605008
## 2009-02-02 0.04082940 0.14484800 -0.1312758 -0.05568857 0.088684084
## d_35 d_36 d_37 d_38 d_39
## 2008-12-29 -0.40281945 -0.38209028 -0.38345643 -0.26410649 -0.32678654
## 2009-01-05 -0.31757678 -0.30158527 -0.30063421 -0.18952899 -0.25676739
## 2009-01-12 -0.15236790 -0.13984462 -0.13732926 -0.03360909 -0.10521652
## 2009-01-19 -0.17758478 -0.17105538 -0.16543749 -0.07040157 -0.14644021
## 2009-01-26 -0.11324000 -0.11102894 -0.10328177 -0.01558557 -0.09564504
## 2009-02-02 -0.02587609 -0.02716283 -0.01773624 0.06355225 -0.02011787
## d_4 d_40 d_41 d_42 d_43
## 2008-12-29 -0.280539857 -0.378186634 -0.36906123 -0.33153862 -0.34417146
## 2009-01-05 -0.193835910 -0.293015403 -0.29328637 -0.24724464 -0.26981802
## 2009-01-12 0.002584089 -0.098100517 -0.13690188 -0.05319120 -0.11411814
## 2009-01-19 -0.059795882 -0.161923794 -0.17164637 -0.11784055 -0.15112161
## 2009-01-26 0.001965910 -0.101509302 -0.11568549 -0.05819719 -0.09650256
## 2009-02-02 0.097571091 -0.007126749 -0.03576155 0.03548559 -0.01754347
## d_44 d_45 d_46 d_47 d_48
## 2008-12-29 -0.31504188 -0.306773132 -0.33081533 -0.34217411 -0.351014741
## 2009-01-05 -0.23591917 -0.236382577 -0.25049440 -0.25879469 -0.270463965
## 2009-01-12 -0.07415809 -0.085285401 -0.08544059 -0.08917624 -0.105184448
## 2009-01-19 -0.10845384 -0.125100094 -0.12134685 -0.12422313 -0.140874272
## 2009-01-26 -0.04982144 -0.073872157 -0.06193603 -0.06232659 -0.081261407
## 2009-02-02 0.03326945 0.001756799 0.02261350 0.02503365 0.003471469
## d_49 d_5 d_50 d_51 d_52
## 2008-12-29 -0.218660197 -0.34347179 -0.31382418 -0.36837408 -0.28965622
## 2009-01-05 -0.149111905 -0.26070363 -0.23805546 -0.28757662 -0.21187081
## 2009-01-12 -0.002829707 -0.09663337 -0.07747187 -0.12673288 -0.05351196
## 2009-01-19 -0.038276033 -0.12585175 -0.11766479 -0.15613565 -0.08636317
## 2009-01-26 0.012729499 -0.06384986 -0.06225548 -0.09569833 -0.02863494
## 2009-02-02 0.086738962 0.02184702 0.01866282 -0.01188030 0.05289279
## d_53 d_54 d_55 d_56 d_57
## 2008-12-29 -0.283092116 -0.32677770 -0.45949961 -0.33793421 -0.39809675
## 2009-01-05 -0.199533249 -0.24646550 -0.37284749 -0.25636081 -0.31974754
## 2009-01-12 -0.006201686 -0.07985904 -0.17647839 -0.09346380 -0.16074970
## 2009-01-19 -0.071543269 -0.11779423 -0.23890717 -0.12380726 -0.19318048
## 2009-01-26 -0.012546088 -0.05859387 -0.17719094 -0.06285560 -0.13496776
## 2009-02-02 0.080550305 0.02631969 -0.08162711 0.02188823 -0.05297001
## d_58 d_6 d_60 d_61 d_62
## 2008-12-29 -0.29602902 -0.341605339 -0.34345968 -0.371973834 -0.37267268
## 2009-01-05 -0.21791059 -0.262766807 -0.26210530 -0.286281982 -0.28595053
## 2009-01-12 -0.05913938 -0.101284785 -0.09447542 -0.120073744 -0.11672687
## 2009-01-19 -0.09178749 -0.135848126 -0.13142922 -0.145590216 -0.14386649
## 2009-01-26 -0.03377763 -0.077465522 -0.07131275 -0.080923111 -0.07855397
## 2009-02-02 0.04803602 0.005398679 0.01443214 0.006931764 0.01059889
## d_63 d_64 d_65 d_66 d_67
## 2008-12-29 -0.37681039 -0.4959411 -0.31895041 -0.33088536 -0.191873712
## 2009-01-05 -0.29518364 -0.4122248 -0.23832811 -0.25300314 -0.116906363
## 2009-01-12 -0.13352556 -0.2480418 -0.07647227 -0.09373081 0.067988421
## 2009-01-19 -0.16214743 -0.2753081 -0.10757258 -0.12755011 -0.005443494
## 2009-01-26 -0.10098115 -0.2123664 -0.04744304 -0.06984309 0.046001510
## 2009-02-02 -0.01650161 -0.1261076 0.03651665 0.01195635 0.132244565
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 -0.36359534 -0.380082266 -0.31851211 NA NA NA NA NA
## 2009-01-05 -0.27708875 -0.295911041 -0.24011260 NA NA NA NA NA
## 2009-01-12 -0.11063872 -0.131196041 -0.08162365 NA NA NA NA NA
## 2009-01-19 -0.13466541 -0.158144434 -0.11328450 NA NA NA NA NA
## 2009-01-26 -0.06920962 -0.094814001 -0.05495505 NA NA NA NA NA
## 2009-02-02 0.01916250 -0.008172109 0.02695015 NA NA NA NA NA
tail(pred.A.cv$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 1.0507262 1.0086899 1.0671375 1.0404536 1.0471231 0.9267457
## 2020-11-09 NA 0.6495543 0.6084700 0.6480466 0.6275711 0.6160250 0.5264845
## 2020-11-16 NA 0.9930112 0.9705669 1.0050999 0.9899644 0.9873549 0.8734524
## 2020-11-23 NA 0.9106816 0.8809822 0.9241371 0.9018572 0.9019613 0.7762845
## 2020-11-30 NA 1.0528323 1.0259840 1.0789589 1.0540887 1.0631642 0.9233131
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2
## 2020-11-02 1.1210214 1.1371417 1.1236092 1.0229126 1.0957606 1.3054840
## 2020-11-09 0.7265924 0.7185737 0.7140727 0.6159464 0.7033632 0.8580272
## 2020-11-16 1.0895570 1.1127847 1.0676951 0.9712165 1.0664872 1.2371962
## 2020-11-23 0.9925457 1.0048637 0.9794156 0.8831923 0.9752908 1.1592842
## 2020-11-30 1.1249627 1.1944760 1.1244620 1.0297807 1.1122199 1.3150344
## 2020-12-07 NA NA NA NA NA NA
## d_20 d_21 d_22 d_23 d_25 d_26
## 2020-11-02 1.0059040 1.0746504 1.1082391 1.0658122 0.9698567 1.0602988
## 2020-11-09 0.6064795 0.6652905 0.6985243 0.6689346 0.5706359 0.6294221
## 2020-11-16 0.9776844 1.0217154 1.0483121 1.0459459 0.9323433 0.9981913
## 2020-11-23 0.8909408 0.9347339 0.9629368 0.9524129 0.8412692 0.9165426
## 2020-11-30 1.0582063 1.0803400 1.1069639 1.0993921 0.9905781 1.0734375
## 2020-12-07 NA NA NA NA NA NA
## d_28 d_29 d_3 d_31 d_32 d_33
## 2020-11-02 0.9248176 0.9534998 1.0844677 0.9818748 0.9623272 0.9882907
## 2020-11-09 0.5331957 0.5551532 0.6645467 0.5863394 0.5569678 0.5951057
## 2020-11-16 0.8833477 0.9043996 1.0318881 0.9542272 0.9117429 0.9368148
## 2020-11-23 0.7952065 0.8149507 0.9422594 0.8663463 0.8209447 0.8550988
## 2020-11-30 0.9398868 0.9579412 1.0953930 1.0129969 0.9667554 0.9939886
## 2020-12-07 NA NA NA NA NA NA
## d_34 d_35 d_36 d_37 d_38 d_39
## 2020-11-02 1.0294080 1.0642775 0.9974888 1.0412300 1.0219983 0.8545244
## 2020-11-09 0.6253148 0.6541375 0.5928055 0.6345602 0.6173857 0.4726909
## 2020-11-16 0.9879691 1.0137144 0.9496459 0.9941513 0.9793151 0.8159157
## 2020-11-23 0.9019640 0.9213251 0.8613282 0.9032006 0.8882420 0.7353018
## 2020-11-30 1.0538837 1.0726021 1.0084492 1.0502736 1.0412946 0.8766090
## 2020-12-07 NA NA NA NA NA NA
## d_4 d_40 d_41 d_42 d_43 d_44
## 2020-11-02 1.0681120 0.7780578 0.9466048 0.9926037 0.9394276 1.0345445
## 2020-11-09 0.6580025 0.4083154 0.5407011 0.5885471 0.5380963 0.6262167
## 2020-11-16 1.0092313 0.7225507 0.9058082 0.9440133 0.8993429 0.9819527
## 2020-11-23 0.9395216 0.6279406 0.8108408 0.8694003 0.8076413 0.8994283
## 2020-11-30 1.0813790 0.7542407 0.9657512 1.0090471 0.9590759 1.0472659
## 2020-12-07 NA NA NA NA NA NA
## d_45 d_46 d_47 d_48 d_49 d_5
## 2020-11-02 0.9096160 1.0072557 1.0539982 1.0031619 1.0166299 1.0012875
## 2020-11-09 0.5115848 0.6123676 0.6645520 0.6063096 0.6102352 0.6105413
## 2020-11-16 0.8716326 0.9556017 1.0225498 0.9557000 0.9820844 0.9428539
## 2020-11-23 0.7825943 0.8774532 0.9356096 0.8730891 0.8911619 0.8586572
## 2020-11-30 0.9329601 1.0189961 1.0618455 1.0152990 1.0380877 0.9989886
## 2020-12-07 NA NA NA NA NA NA
## d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 0.9410239 1.0137597 1.0508521 1.1625364 1.0149022 0.8914866
## 2020-11-09 0.5581723 0.6196333 0.6393464 0.7178098 0.6139069 0.4989899
## 2020-11-16 0.8968702 0.9809427 0.9984167 1.0876010 0.9784398 0.8483559
## 2020-11-23 0.8192968 0.8897444 0.9143677 1.0129602 0.8936533 0.7700250
## 2020-11-30 0.9576246 1.0240487 1.0644770 1.1564030 1.0519860 0.9058540
## 2020-12-07 NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61
## 2020-11-02 1.0128975 0.9317901 1.0767284 1.0112969 1.0041229 1.0596957
## 2020-11-09 0.6145393 0.5403825 0.6649401 0.6046561 0.6029944 0.6521355
## 2020-11-16 0.9535509 0.8878284 1.0268115 0.9653383 0.9527762 0.9958379
## 2020-11-23 0.8740066 0.8018455 0.9409341 0.8776155 0.8655090 0.9144330
## 2020-11-30 1.0179338 0.9448624 1.0858201 1.0249889 1.0113756 1.0610470
## 2020-12-07 NA NA NA NA NA NA
## d_62 d_63 d_64 d_65 d_66 d_67
## 2020-11-02 1.0578885 1.0221486 0.9190848 1.0655690 1.0020878 1.0011043
## 2020-11-09 0.6563047 0.6177055 0.5248825 0.6279024 0.5975990 0.6023341
## 2020-11-16 0.9949120 0.9714508 0.8706186 0.9945009 0.9557804 0.9743314
## 2020-11-23 0.9145625 0.8853176 0.7833423 0.9168700 0.8659570 0.8928677
## 2020-11-30 1.0577934 1.0312216 0.9255889 1.0788609 1.0167814 1.0280331
## 2020-12-07 NA NA NA NA NA NA
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.1669645 1.0304772 1.0327063 NA NA NA NA NA
## 2020-11-09 0.7444754 0.6205897 0.6260371 NA NA NA NA NA
## 2020-11-16 1.1049439 0.9682406 0.9829630 NA NA NA NA NA
## 2020-11-23 1.0173789 0.8835628 0.8977986 NA NA NA NA NA
## 2020-11-30 1.1628995 1.0319215 1.0454443 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA NA NA
head(pred.A.cv.log$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 0.6706146 0.7606878 0.7965575 0.7612359 0.7728678 0.6447413
## 2009-01-05 NA 0.7318498 0.8202008 0.8606169 0.8214333 0.8325421 0.7002994
## 2009-01-12 NA 0.8662736 0.9587303 1.0078974 0.9570165 0.9683108 0.8211545
## 2009-01-19 NA 0.8450322 0.9255531 0.9750002 0.9294927 0.9389298 0.8024734
## 2009-01-26 NA 0.9027973 0.9785070 1.0326554 0.9839063 0.9923840 0.8543571
## 2009-02-02 NA 0.9872792 1.0596663 1.1201223 1.0651955 1.0728906 0.9297896
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2008-12-29 0.6946862 1.259380 0.6883669 0.6813451 0.7848421 1.309593 0.8131327
## 2009-01-05 0.7600640 1.320701 0.7522142 0.7406528 0.8510093 1.378013 0.8748953
## 2009-01-12 0.9049800 1.496132 0.8903847 0.8732426 1.0059346 1.575808 1.0204620
## 2009-01-19 0.8810218 1.414538 0.8710917 0.8488276 0.9691096 1.482930 0.9833135
## 2009-01-26 0.9429111 1.460218 0.9318822 0.9037710 1.0288504 1.534191 1.0376768
## 2009-02-02 1.0340324 1.545246 1.0198560 0.9852654 1.1206796 1.630296 1.1218578
## d_21 d_22 d_23 d_25 d_26 d_28
## 2008-12-29 0.7005308 0.6641760 0.9899413 0.7424012 0.7526664 0.6861555
## 2009-01-05 0.7622531 0.7273895 1.0529158 0.8007109 0.8140902 0.7424850
## 2009-01-12 0.8990056 0.8651092 1.2094413 0.9400527 0.9564450 0.8745136
## 2009-01-19 0.8753760 0.8453610 1.1587992 0.9028977 0.9237048 0.8425539
## 2009-01-26 0.9329034 0.9058720 1.2111710 0.9542259 0.9790364 0.8930292
## 2009-02-02 1.0176520 0.9938189 1.2962252 1.0345452 1.0634588 0.9707403
## d_29 d_3 d_31 d_32 d_33 d_34
## 2008-12-29 0.6124998 0.7551884 0.8974433 0.6137073 0.6641565 0.8094408
## 2009-01-05 0.6676277 0.8168745 0.9564308 0.6689308 0.7235239 0.8716082
## 2009-01-12 0.7887518 0.9539608 1.1060374 0.7907745 0.8591761 1.0232820
## 2009-01-19 0.7692809 0.9286305 1.0550829 0.7706217 0.8313757 0.9793967
## 2009-01-26 0.8210921 0.9850788 1.1039452 0.8224522 0.8862606 1.0334441
## 2009-02-02 0.8969287 1.0685177 1.1846203 0.8985467 0.9689700 1.1194652
## d_35 d_36 d_37 d_38 d_39 d_4
## 2008-12-29 0.6857350 0.7001492 0.6996790 0.7878260 0.7398876 0.7752066
## 2009-01-05 0.7465360 0.7585759 0.7597895 0.8485247 0.7933479 0.8451613
## 2009-01-12 0.8804185 0.8914687 0.8942559 0.9913879 0.9229616 1.0283162
## 2009-01-19 0.8583079 0.8638445 0.8692075 0.9553203 0.8855175 0.9659078
## 2009-01-26 0.9151873 0.9170848 0.9247184 1.0089276 0.9315126 1.0272508
## 2009-02-02 0.9985963 0.9971405 1.0071074 1.0918260 1.0044663 1.1301398
## d_40 d_41 d_42 d_43 d_44 d_45
## 2008-12-29 0.7030886 0.7098239 0.7366633 0.7272078 0.7502807 0.7554435
## 2009-01-05 0.7653613 0.7653929 0.8012066 0.7830606 0.8116354 0.8102097
## 2009-01-12 0.9298222 0.8946382 0.9725317 0.9147006 0.9537060 0.9420288
## 2009-01-19 0.8721318 0.8638274 0.9114382 0.8812370 0.9211939 0.9049858
## 2009-01-26 0.9262704 0.9133193 0.9672706 0.9305039 0.9765080 0.9523179
## 2009-02-02 1.0178001 0.9891166 1.0621082 1.0067787 1.0608435 1.0269316
## d_46 d_47 d_48 d_49 d_5 d_50
## 2008-12-29 0.7360538 0.7282917 0.7213351 0.8257199 0.7276448 0.7486671
## 2009-01-05 0.7975489 0.7914178 0.7817802 0.8847847 0.7902313 0.8075315
## 2009-01-12 0.9406072 0.9375028 0.9222181 1.0237441 0.9309176 0.9481325
## 2009-01-19 0.9073775 0.9050421 0.8898306 0.9877512 0.9039379 0.9107246
## 2009-01-26 0.9628709 0.9626796 0.9444417 1.0391444 0.9616061 0.9625634
## 2009-02-02 1.0477803 1.0504310 1.0279144 1.1187164 1.0475156 1.0436490
## d_51 d_52 d_53 d_54 d_55 d_56
## 2008-12-29 0.7097669 0.7684857 0.7732306 0.7395916 0.6481809 0.7316854
## 2009-01-05 0.7692716 0.8303148 0.8403598 0.8012357 0.7066362 0.7936705
## 2009-01-12 0.9032798 0.9724409 1.0193212 0.9462785 0.8597275 0.9338729
## 2009-01-19 0.8769167 0.9407300 0.9546273 0.9108793 0.8075114 0.9057879
## 2009-01-26 0.9313829 0.9963874 1.0124510 0.9662798 0.8587558 0.9625627
## 2009-02-02 1.0126707 1.0808106 1.1110667 1.0517828 0.9447293 1.0475588
## d_57 d_58 d_6 d_60 d_61 d_62
## 2008-12-29 0.6890545 0.7630993 0.7306130 0.7273561 0.7072918 0.7082638
## 2009-01-05 0.7449823 0.8248466 0.7901349 0.7888020 0.7703354 0.7720273
## 2009-01-12 0.8731314 0.9665146 0.9281829 0.9325480 0.9093790 0.9139599
## 2009-01-19 0.8450733 0.9352523 0.8963010 0.8985438 0.8862629 0.8891429
## 2009-01-26 0.8955566 0.9909217 0.9498831 0.9540676 0.9452884 0.9488497
## 2009-02-02 0.9719376 1.0752384 1.0316852 1.0393537 1.0319375 1.0370641
## d_63 d_64 d_65 d_66 d_67 d_68
## 2008-12-29 0.7038043 0.6252409 0.7457820 0.7368612 0.8470805 0.7131668
## 2009-01-05 0.7634419 0.6795636 0.8081115 0.7963399 0.9127455 0.7773820
## 2009-01-12 0.8971649 0.8005343 0.9497917 0.9336236 1.0978207 0.9179349
## 2009-01-19 0.8716607 0.7787665 0.9204619 0.9024040 1.0198602 0.8959479
## 2009-01-26 0.9264755 0.8291530 0.9772923 0.9558602 1.0734972 0.9563837
## 2009-02-02 1.0080016 0.9036717 1.0627033 1.0372061 1.1700129 1.0445998
## d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.7015800 0.7460545 NA NA NA NA NA
## 2009-01-05 0.7629534 0.8066647 NA NA NA NA NA
## 2009-01-12 0.8993207 0.9449591 NA NA NA NA NA
## 2009-01-19 0.8752061 0.9153104 NA NA NA NA NA
## 2009-01-26 0.9322483 0.9701141 NA NA NA NA NA
## 2009-02-02 1.0164684 1.0527666 NA NA NA NA NA
tail(pred.A.cv.log$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 2.882844 2.758274 2.927953 2.849952 2.866982 2.543824
## 2020-11-09 NA 1.930042 1.848378 1.925364 1.885719 1.862744 1.704543
## 2020-11-16 NA 2.720860 2.654734 2.751349 2.709094 2.700110 2.411326
## 2020-11-23 NA 2.505726 2.427134 2.537214 2.480444 2.478935 2.187893
## 2020-11-30 NA 2.888383 2.805759 2.961924 2.888140 2.912394 2.534285
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 3.082275 3.132966 3.096215 2.798747 3.006990 3.711721 2.749532
## 2020-11-09 2.077604 2.061225 2.055608 1.862871 2.030897 2.372429 1.843946
## 2020-11-16 2.986676 3.056966 2.927466 2.657325 2.919903 3.465943 2.672555
## 2020-11-23 2.710504 2.744040 2.679987 2.433282 2.665290 3.205905 2.450346
## 2020-11-30 3.094232 3.316756 3.098202 2.817318 3.056313 3.745974 2.896332
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2020-11-02 2.948014 3.055046 2.918045 2.652956 2.907973 2.538925 2.613302
## 2020-11-09 1.957553 2.027810 1.961922 1.779674 1.889771 1.716171 1.754513
## 2020-11-16 2.795636 2.876708 2.860074 2.555170 2.732254 2.435688 2.487747
## 2020-11-23 2.562619 2.641095 2.604509 2.332708 2.517836 2.230161 2.274773
## 2020-11-30 2.964170 3.050051 3.016710 2.708316 2.945363 2.577299 2.624347
## 2020-12-07 NA NA NA NA NA NA NA
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2020-11-02 2.976240 2.684562 2.633946 2.708008 2.815652 2.915474 2.728219
## 2020-11-09 1.955466 1.807451 1.755997 1.827529 1.879557 1.934444 1.820077
## 2020-11-16 2.823231 2.611054 2.503633 2.571844 2.701046 2.771363 2.600361
## 2020-11-23 2.581016 2.391288 2.286197 2.369946 2.478351 2.526675 2.380422
## 2020-11-30 3.007955 2.768900 2.644962 2.722974 2.884871 2.939221 2.757583
## 2020-12-07 NA NA NA NA NA NA NA
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2020-11-02 2.850910 2.794233 2.370156 2.934581 2.193477 2.592523 2.719307
## 2020-11-09 1.898133 1.864248 1.617786 1.947176 1.515386 1.727421 1.815290
## 2020-11-16 2.719326 2.677057 2.280132 2.766407 2.074761 2.488447 2.589984
## 2020-11-23 2.482762 2.443888 2.103451 2.580007 1.887377 2.262860 2.403657
## 2020-11-30 2.875985 2.847947 2.422631 2.973119 2.141378 2.641879 2.763781
## 2020-12-07 NA NA NA NA NA NA NA
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2020-11-02 2.573452 2.837354 2.500204 2.760023 2.884754 2.746439 2.779292
## 2020-11-09 1.722592 1.885929 1.679076 1.859538 1.954098 1.846754 1.850934
## 2020-11-16 2.471951 2.691393 2.406597 2.620969 2.795120 2.619024 2.684384
## 2020-11-23 2.255228 2.478012 2.201445 2.423908 2.562271 2.411325 2.450908
## 2020-11-30 2.623847 2.872641 2.558523 2.792431 2.906929 2.779788 2.838648
## 2020-12-07 NA NA NA NA NA NA NA
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 2.744472 2.584436 2.770826 2.879907 3.216428 2.774320 2.458789
## 2020-11-09 1.856660 1.762323 1.868146 1.908189 2.061572 1.857712 1.660464
## 2020-11-16 2.588403 2.472705 2.681023 2.732310 2.983805 2.674674 2.354676
## 2020-11-23 2.379295 2.288108 2.447224 2.511893 2.769067 2.457145 2.177168
## 2020-11-30 2.737660 2.627523 2.798887 2.918579 3.196047 2.878588 2.493821
## 2020-12-07 NA NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 2.776892 2.558377 2.952317 2.769914 2.747011 2.908258 2.908396
## 2020-11-09 1.864346 1.729599 1.955655 1.844211 1.839182 1.934634 1.946232
## 2020-11-16 2.616589 2.448009 2.808177 2.644907 2.609221 2.727976 2.730283
## 2020-11-23 2.416420 2.246209 2.576958 2.422584 2.391077 2.514584 2.519293
## 2020-11-30 2.790393 2.591462 2.978609 2.807082 2.766486 2.911542 2.907074
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 2.798118 2.526940 2.920288 2.740792 2.736374 3.230491 2.822425
## 2020-11-09 1.867185 1.703542 1.885002 1.828861 1.836364 2.117154 1.873171
## 2020-11-16 2.659452 2.406970 2.719528 2.616464 2.663721 3.035826 2.651759
## 2020-11-23 2.439862 2.205669 2.516257 2.391593 2.455210 2.781173 2.436340
## 2020-11-30 2.823025 2.542705 2.958610 2.780836 2.810434 3.216702 2.825872
## 2020-12-07 NA NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.827221 NA NA NA NA NA
## 2020-11-09 1.882410 NA NA NA NA NA
## 2020-11-16 2.689678 NA NA NA NA NA
## 2020-11-23 2.469984 NA NA NA NA NA
## 2020-11-30 2.862860 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA
summary(pred.A.cv)
## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 712 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX
## obs 0.16720556 0.16720556 0.11062431
## average 0.09258214 0.09258214 0.04833502
##
## R2:
## EX.mu EX.mu.beta EX
## obs 0.5990090 0.5990090 0.8244769
## average 0.6215252 0.6215252 0.8968412
##
## Coverage of 95% prediction intervals:
## EX
## obs 0.9536517
## average 0.9047619
summary(pred.A.cv.log)
## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 712 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.2211963 0.2211963 0.15552805 0.15542748
## average 0.1045434 0.1045434 0.05577006 0.05597244
##
## R2:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.6185378 0.6185378 0.8114124 0.8116562
## average 0.6822837 0.6822837 0.9095833 0.9089259
##
## Coverage of 95% prediction intervals:
## EX.pred
## obs 0.9536517
## average 0.9523810
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModA.jpeg"),
width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()
## quartz_off_screen
## 2
What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites
x.A <- coef(est.denver.model.A, pars = "cov")$par
x.A
## [1] -7.558680 -7.132387 12.520027 -3.032854 -4.759856
E.A <- predict(denver.model.A, est.denver.model.A, LTA = T,
transform="unbiased", pred.var=FALSE)
print(E.A)
## Prediction for STmodel.
##
## Regression parameters:
## 0 Spatio-temporal covariate(s).
## 16 beta-fields regression parameters in x$pars.
##
## Regression parameters are assumed to be known and
## prediction variances do NOT include
## uncertainties in regression parameters.
##
## Prediction of beta-fields, (x$beta):
## List of 2
## $ mu: num [1:69, 1:2] 0.2211736 0.0025188 -0.0001717 0.0546205 0.0000313 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX: num [1:69, 1:2] 0.23143 0.01497 -0.00252 0.03873 -0.00599 ...
## ..- attr(*, "dimnames")=List of 2
##
## Predictions for log-Gaussian field of type: unbiased
##
## Predictions for 624 times at 69 locations.
## List of 5
## $ EX.mu : num [1:624, 1:69] 1.19 1.25 1.43 1.35 1.39 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.mu.beta: num [1:624, 1:69] 1.25 1.31 1.49 1.4 1.44 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX : num [1:624, 1:69] 1.28 1.34 1.53 1.43 1.48 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.pred : num [1:624, 1:69] 1.29 1.35 1.53 1.44 1.48 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.EX : num [1:624, 1:69] 0.223 0.27 0.399 0.336 0.367 ...
## ..- attr(*, "dimnames")=List of 2
##
## Variances have been computed.
## List of 2
## $ log.VX : num [1:624, 1:69] 0.0487 0.0486 0.0486 0.0485 0.0485 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.VX.pred: num [1:624, 1:69] 0.0573 0.0572 0.0571 0.0571 0.057 ...
## ..- attr(*, "dimnames")=List of 2
##
## Mean squared prediciton errors NOT been computed.
##
## 69 temporal averages have been compute.
## List of 1
## $ LTA:'data.frame': 69 obs. of 4 variables:
week_preds <- as.data.frame(E.A$EX) %>%
mutate(week = as.Date(rownames(E.A$EX)),
year = year(as.Date(rownames(E.A$EX))))
week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]
box_preds <- week_preds %>%
select(week, all_of(week_sites)) %>%
#filter(week %in% week_weeks) %>%
mutate(month = month(week),
year = year(week)) %>%
filter(year %in% c(2009:2020))
box_data <- box_preds %>%
pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)
## week month year location
## Min. :2009-01-05 Min. : 1.000 Min. :2009 Length:42364
## 1st Qu.:2011-12-26 1st Qu.: 4.000 1st Qu.:2011 Class :character
## Median :2014-12-22 Median : 7.000 Median :2014 Mode :character
## Mean :2014-12-22 Mean : 6.506 Mean :2014
## 3rd Qu.:2017-12-18 3rd Qu.: 9.000 3rd Qu.:2017
## Max. :2020-12-07 Max. :12.000 Max. :2020
##
## pred
## Min. :0.3344
## 1st Qu.:1.1149
## Median :1.3724
## Mean :1.4213
## 3rd Qu.:1.6580
## Max. :3.8895
## NA's :1080
#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
group_by(month) %>%
summarize(median = median(pred, na.rm=T)) %>%
arrange(desc(median))
box_summary
pred_box_plot <- ggplot(box_data) +
geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
show.legend = F) +
#scale_color_viridis(name = "Month", discrete = T) +
facet_wrap(. ~ year, scales = "free") +
xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
# theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot
For this version of the model, use iid for cov.beta (beta, beta1, beta2) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.
denver.data.B <- denver.data2.2
names(denver.data.B$covars)
## [1] "ID" "lon" "lat" "impervious_2500"
## [5] "open_2500" "med_int_2500" "high_int_100" "dist_m_cafos"
## [9] "aadt_100" "aadt_250" "x" "y"
## [13] "type"
LUR.B <- list(covar_fun, covar_fun, covar_fun)
cov.beta.B <- list(covf = c("iid", "iid", "iid"), nugget = T)
cov.nu.B <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.B <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")
denver.model.B <- createSTmodel(denver.data.B, LUR = LUR.B,
ST = c("bc_st_no2", "bc_st_smk"),
cov.beta = cov.beta.B, cov.nu = cov.nu.B,
locations = locations.B)
## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.
denver.model.B
## STmodel-object with:
## No. locations: 69 (observed: 64)
## No. time points: 624 (observed: 231)
## No. obs: 943
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## Models for the beta-fields are:
## $const
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 +
## dist_m_cafos + aadt_100 + aadt_250
##
## $V1
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 +
## dist_m_cafos + aadt_100 + aadt_250
##
## $V2
## ~~impervious_2500 + open_2500 + med_int_2500 + high_int_100 +
## dist_m_cafos + aadt_100 + aadt_250
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## Covariance model for the beta-field(s):
## Covariance type(s): iid, iid, iid
## Nugget: Yes, Yes, Yes
## Covariance model for the nu-field(s):
## Covariance type: exp
## Nugget: ~1
## Random effect: No
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 712
## Dates: 2018-05-07 to 2020-04-06
names <- loglikeSTnames(denver.model.B, all=FALSE)
names
## [1] "log.nugget.const.iid" "log.nugget.V1.iid"
## [3] "log.nugget.V2.iid" "nu.log.range.exp"
## [5] "nu.log.sill.exp" "nu.log.nugget.(Intercept).exp"
x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
c(-10, -5, -5, 8, -3, -5),
c(-10, -5, -5, 8, -5, -5),
c(-10, -5, -5, 10, -3, -5),
c(-10, -5, -5, 10, -5, -5),
c(-12, -5, -5, 8, -5, -5),
c(-12, -5, -5, 10, -5, -5),
c(-12, -5, -5, 12, -5, -5))
rownames(x.init.B) <- loglikeSTnames(denver.model.B, all=FALSE)
x.init.B
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid 0 -10 -10 -10 -10 -12 -12 -12
## log.nugget.V1.iid 0 -5 -5 -5 -5 -5 -5 -5
## log.nugget.V2.iid 0 -5 -5 -5 -5 -5 -5 -5
## nu.log.range.exp 0 8 8 10 10 8 10 12
## nu.log.sill.exp 0 -3 -5 -3 -5 -5 -5 -5
## nu.log.nugget.(Intercept).exp 0 -5 -5 -5 -5 -5 -5 -5
#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.B <- estimate.STmodel(denver.model.B, x.init.B)
## Optimisation using starting value 1/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 474.52 |proj g|= 15
## At iterate 10 f = -1257.9 |proj g|= 0.62776
## At iterate 20 f = -1258.2 |proj g|= 0.34897
##
## iterations 29
## function evaluations 35
## segments explored during Cauchy searches 34
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.000608333
## final function value -1258.28
##
## F = -1258.28
## final value -1258.282191
## converged
## Optimisation using starting value 2/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1145.5 |proj g|= 12
## At iterate 10 f = -1458.9 |proj g|= 5.2642
## At iterate 20 f = -1459.4 |proj g|= 0.069819
## At iterate 30 f = -1459.5 |proj g|= 0.85468
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1459.490471
## stopped after 39 iterations
## Optimisation using starting value 3/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1143.1 |proj g|= 20
## At iterate 10 f = -1456.1 |proj g|= 3.879
## At iterate 20 f = -1459.4 |proj g|= 0.14248
## At iterate 30 f = -1459.4 |proj g|= 1.1349
## At iterate 40 f = -1459.5 |proj g|= 0.072186
##
## iterations 48
## function evaluations 61
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0304898
## final function value -1459.49
##
## F = -1459.49
## final value -1459.490601
## converged
## Optimisation using starting value 4/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1309.4 |proj g|= 12
## At iterate 10 f = -1458.9 |proj g|= 11.518
## At iterate 20 f = -1459.4 |proj g|= 0.96479
## At iterate 30 f = -1459.5 |proj g|= 0.031888
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1459.490603
## stopped after 30 iterations
## Optimisation using starting value 5/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1200.7 |proj g|= 20
## At iterate 10 f = -1459.2 |proj g|= 1.2752
## At iterate 20 f = -1459.4 |proj g|= 0.078369
## At iterate 30 f = -1459.4 |proj g|= 0.56669
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1459.490611
## stopped after 38 iterations
## Optimisation using starting value 6/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1142.8 |proj g|= 20
## At iterate 10 f = -1456.6 |proj g|= 10.275
## At iterate 20 f = -1459.3 |proj g|= 0.037973
##
## iterations 25
## function evaluations 40
## segments explored during Cauchy searches 30
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0493346
## final function value -1459.35
##
## F = -1459.35
## final value -1459.348069
## converged
## Optimisation using starting value 7/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1200.4 |proj g|= 20
## At iterate 10 f = -1459.2 |proj g|= 1.2461
## At iterate 20 f = -1459.3 |proj g|= 0.028556
##
## iterations 20
## function evaluations 36
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0285564
## final function value -1459.35
##
## F = -1459.35
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1459.347897
## converged
## Optimisation using starting value 8/8
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1200 |proj g|= 20
## At iterate 10 f = -1459 |proj g|= 4.961
## At iterate 20 f = -1459.3 |proj g|= 0.026597
## At iterate 30 f = -1459.3 |proj g|= 0.4351
## ys=-7.527e-03 -gs= 1.809e-02, BFGS update SKIPPED
## At iterate 40 f = -1459.5 |proj g|= 0.24613
##
## iterations 46
## function evaluations 63
## segments explored during Cauchy searches 50
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0258296
## final function value -1459.49
##
## F = -1459.49
## final value -1459.490552
## converged
print(est.denver.model.B)
## Optimisation for STmodel with 8 starting points.
## Results: 2 converged, 6 not converged, 0 failed.
## Best result for starting point 5, optimisation has NOT converged
## Best converged result for starting point 3
##
## No fixed parameters.
##
## Estimated parameters for all starting point(s):
## [,1] [,2] [,3]
## gamma.bc_st_no2 0.006825267 0.0086902950 0.0086810932
## gamma.bc_st_smk 0.050977166 0.0490210509 0.0490347152
## alpha.const.(Intercept) 0.175407462 0.1956111940 0.1958073526
## alpha.const.impervious_2500 0.002096213 0.0152705785 0.0153226189
## alpha.const.open_2500 0.007918598 0.0064206957 0.0064102731
## alpha.const.med_int_2500 0.016339133 -0.0003250761 -0.0003660604
## alpha.const.high_int_100 0.002482419 0.0048685312 0.0048472442
## alpha.const.dist_m_cafos -0.021717818 -0.0214981651 -0.0215274974
## alpha.const.aadt_100 0.030277467 0.0166623609 0.0166407242
## alpha.const.aadt_250 0.022896805 0.0178550584 0.0178817712
## alpha.V1.(Intercept) 0.166162945 0.1783673742 0.1783837079
## alpha.V1.impervious_2500 0.014938297 0.0075446851 0.0075584878
## alpha.V1.open_2500 -0.006466654 0.0006784633 0.0006804785
## alpha.V1.med_int_2500 -0.007280071 -0.0036866182 -0.0036926297
## alpha.V1.high_int_100 0.023599868 0.0238944621 0.0238902363
## alpha.V1.dist_m_cafos 0.004421040 -0.0024551014 -0.0024314856
## alpha.V1.aadt_100 -0.015712578 -0.0206936558 -0.0206953131
## alpha.V1.aadt_250 -0.031984221 -0.0200627761 -0.0200603014
## alpha.V2.(Intercept) 0.157566420 0.1723380449 0.1723940701
## alpha.V2.impervious_2500 -0.022409303 -0.0170009993 -0.0169686559
## alpha.V2.open_2500 0.013478857 0.0117515720 0.0117436817
## alpha.V2.med_int_2500 0.026519369 0.0178107329 0.0177856997
## alpha.V2.high_int_100 -0.006151440 -0.0041084710 -0.0041163634
## alpha.V2.dist_m_cafos 0.010196517 0.0100957490 0.0100791610
## alpha.V2.aadt_100 -0.021198667 -0.0268102887 -0.0268312866
## alpha.V2.aadt_250 0.004345172 0.0021175995 0.0021302273
## log.nugget.const.iid -14.319024162 -8.7133077165 -8.6933320465
## log.nugget.V1.iid -15.000000000 -7.0772541944 -7.0746701234
## log.nugget.V2.iid -15.000000000 -7.5629631080 -7.5683449287
## nu.log.range.exp 0.000000000 13.8028190119 13.8120753074
## nu.log.sill.exp -4.623411406 -3.0088155216 -3.0085177177
## nu.log.nugget.(Intercept).exp -4.149419586 -4.8187117583 -4.8182150293
## [,4] [,5] [,6]
## gamma.bc_st_no2 0.0086853418 0.0086804968 0.0089507205
## gamma.bc_st_smk 0.0490259149 0.0490321152 0.0490594641
## alpha.const.(Intercept) 0.1956978301 0.1958058811 0.1922766174
## alpha.const.impervious_2500 0.0153057073 0.0153252484 0.0140527012
## alpha.const.open_2500 0.0064092579 0.0064068806 0.0069662685
## alpha.const.med_int_2500 -0.0003602946 -0.0003741714 0.0012286750
## alpha.const.high_int_100 0.0048541069 0.0048462292 0.0053625754
## alpha.const.dist_m_cafos -0.0215188945 -0.0215273574 -0.0209388328
## alpha.const.aadt_100 0.0166421541 0.0166374468 0.0174753357
## alpha.const.aadt_250 0.0178793006 0.0178854015 0.0169291157
## alpha.V1.(Intercept) 0.1783822229 0.1783840296 0.1777210937
## alpha.V1.impervious_2500 0.0075506172 0.0075523751 0.0077571785
## alpha.V1.open_2500 0.0006825803 0.0006810681 0.0005800503
## alpha.V1.med_int_2500 -0.0036849352 -0.0036871961 -0.0040671677
## alpha.V1.high_int_100 0.0238863602 0.0238876258 0.0242013537
## alpha.V1.dist_m_cafos -0.0024432810 -0.0024368070 -0.0024536176
## alpha.V1.aadt_100 -0.0206952794 -0.0206947608 -0.0206231665
## alpha.V1.aadt_250 -0.0200593952 -0.0200573132 -0.0203550796
## alpha.V2.(Intercept) 0.1723653835 0.1723915816 0.1717389221
## alpha.V2.impervious_2500 -0.0169806291 -0.0169680165 -0.0176446886
## alpha.V2.open_2500 0.0117443581 0.0117413398 0.0121515502
## alpha.V2.med_int_2500 0.0177917264 0.0177816580 0.0186519515
## alpha.V2.high_int_100 -0.0041138614 -0.0041172903 -0.0038791942
## alpha.V2.dist_m_cafos 0.0100859408 0.0100788070 0.0104441570
## alpha.V2.aadt_100 -0.0268276554 -0.0268338080 -0.0261199038
## alpha.V2.aadt_250 0.0021296855 0.0021331199 0.0015899649
## log.nugget.const.iid -8.6940671217 -8.6870154341 -11.9673583758
## log.nugget.V1.iid -7.0766936210 -7.0748494524 -7.0723990876
## log.nugget.V2.iid -7.5678727748 -7.5681824916 -7.4631928925
## nu.log.range.exp 13.8060928682 13.8095411299 13.8321924769
## nu.log.sill.exp -3.0090198820 -3.0084829692 -3.0118696588
## nu.log.nugget.(Intercept).exp -4.8184536221 -4.8183782834 -4.8080020534
## [,7] [,8]
## gamma.bc_st_no2 0.0089512000 0.0086805587
## gamma.bc_st_smk 0.0490578896 0.0490258663
## alpha.const.(Intercept) 0.1922668708 0.1957613345
## alpha.const.impervious_2500 0.0140465974 0.0153274410
## alpha.const.open_2500 0.0069623102 0.0063987275
## alpha.const.med_int_2500 0.0012280383 -0.0003880943
## alpha.const.high_int_100 0.0053654093 0.0048444566
## alpha.const.dist_m_cafos -0.0209381425 -0.0215300089
## alpha.const.aadt_100 0.0174719150 0.0166265606
## alpha.const.aadt_250 0.0169329549 0.0179017693
## alpha.V1.(Intercept) 0.1777247001 0.1784124588
## alpha.V1.impervious_2500 0.0077603418 0.0075465152
## alpha.V1.open_2500 0.0005855862 0.0006816656
## alpha.V1.med_int_2500 -0.0040655491 -0.0036764101
## alpha.V1.high_int_100 0.0241961668 0.0238800331
## alpha.V1.dist_m_cafos -0.0024543853 -0.0024389340
## alpha.V1.aadt_100 -0.0206214482 -0.0207037662
## alpha.V1.aadt_250 -0.0203575991 -0.0200540593
## alpha.V2.(Intercept) 0.1717365101 0.1723802107
## alpha.V2.impervious_2500 -0.0176526009 -0.0169684392
## alpha.V2.open_2500 0.0121509475 0.0117355998
## alpha.V2.med_int_2500 0.0186569176 0.0177763912
## alpha.V2.high_int_100 -0.0038770257 -0.0041188040
## alpha.V2.dist_m_cafos 0.0104490010 0.0100798319
## alpha.V2.aadt_100 -0.0261205759 -0.0268426034
## alpha.V2.aadt_250 0.0015905254 0.0021414966
## log.nugget.const.iid -11.9786223678 -8.6744466377
## log.nugget.V1.iid -7.0738093944 -7.0784616664
## log.nugget.V2.iid -7.4661100356 -7.5700023965
## nu.log.range.exp 13.8307608958 13.8069718049
## nu.log.sill.exp -3.0117132930 -3.0090709173
## nu.log.nugget.(Intercept).exp -4.8080881154 -4.8185034029
##
## Function value(s):
## [1] 1258.282 1459.490 1459.491 1459.491 1459.491 1459.348 1459.348 1459.491
Define the CV groups
set.seed(1000)
site_idsB <- unique(denver.model.B$obs$ID[which(str_detect(denver.model.B$obs$ID, "d_") == T)])
site_idsB
## [1] "d_1" "d_18" "d_22" "d_25" "d_28" "d_31" "d_33" "d_36" "d_39" "d_4"
## [11] "d_7" "d_11" "d_12" "d_19" "d_23" "d_26" "d_29" "d_32" "d_35" "d_38"
## [21] "d_40" "d_8" "d_42" "d_44" "d_46" "d_50" "d_52" "d_54" "d_56" "d_58"
## [31] "d_61" "d_63" "d_65" "d_68" "d_15" "d_41" "d_43" "d_45" "d_47" "d_51"
## [41] "d_53" "d_67" "d_49" "d_14" "d_48" "d_37" "d_55" "d_57" "d_64" "d_66"
## [51] "d_16" "d_2" "d_5" "d_60" "d_62" "d_34" "d_10" "d_13" "d_21" "d_3"
## [61] "d_6" "d_17" "d_20"
Ind.cv.B <- createCV(denver.model.B, groups = 10, #min.dist = .1,
subset = site_idsB)
ID.cv.B <- sapply(split(denver.model.B$obs$ID, Ind.cv.B), unique)
print(sapply(ID.cv.B, length))
## 0 1 2 3 4 5 6 7 8 9 10
## 1 7 7 7 6 6 6 6 6 6 6
table(Ind.cv.B)
## Ind.cv.B
## 0 1 2 3 4 5 6 7 8 9 10
## 231 66 76 82 67 57 85 72 78 54 75
I.col.B <- apply(sapply(ID.cv.B, function(x) denver.model.B$locations$ID%in% x), 1,
function(x) if(sum(x)==1) which(x) else 0)
names(I.col.B) <- denver.model.B$locations$ID
print(I.col.B)
## central d_1 d_10 d_11 d_12 d_13 d_14 d_15 d_16 d_17
## 1 6 2 3 4 4 4 9 4 10
## d_18 d_19 d_2 d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 8 5 11 3 2 11 4 9 11 9
## d_29 d_3 d_31 d_32 d_33 d_34 d_35 d_36 d_37 d_38
## 2 4 6 8 5 5 10 8 3 8
## d_39 d_4 d_40 d_41 d_42 d_43 d_44 d_45 d_46 d_47
## 6 7 7 3 7 8 11 3 9 5
## d_48 d_49 d_5 d_50 d_51 d_52 d_53 d_54 d_55 d_56
## 9 4 6 9 10 3 7 5 7 6
## d_57 d_58 d_6 d_60 d_61 d_62 d_63 d_64 d_65 d_66
## 2 2 11 5 2 11 10 3 8 6
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 7 10 2 10 0 0 0 0 0
par(mfrow=c(1,1))
plot(denver.model.B$locations$long,
denver.model.B$locations$lat,
pch=23+floor(I.col.B/max(I.col.B)+.5), bg=I.col.B,
xlab="Longitude", ylab="Latitude")
map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)
## [[1]]
## NULL
ID starting conditions using previous model without CV:
x.init.B.cv <- coef(est.denver.model.B, pars="cov")[,c("par","init")]
x.init.B.cv
Run the model with cross validation
est.denver.B.cv <- estimateCV(denver.model.B, x.init.B.cv, Ind.cv.B)
##
## ***************************
## Estimation of CV-set 1/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1325.1 |proj g|= 14.284
## At iterate 10 f = -1325.5 |proj g|= 0.089734
## At iterate 20 f = -1325.6 |proj g|= 0.071136
##
## iterations 25
## function evaluations 41
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0380497
## final function value -1325.55
##
## F = -1325.55
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1325.550743
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1020.4 |proj g|= 20
## At iterate 10 f = -1323.1 |proj g|= 3.4163
## At iterate 20 f = -1325.5 |proj g|= 0.22009
## At iterate 30 f = -1325.5 |proj g|= 0.21144
## At iterate 40 f = -1325.6 |proj g|= 0.028442
##
## iterations 41
## function evaluations 50
## segments explored during Cauchy searches 46
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0159431
## final function value -1325.55
##
## F = -1325.55
## final value -1325.550741
## converged
##
##
## ***************************
## Estimation of CV-set 2/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1337.4 |proj g|= 10.182
## At iterate 10 f = -1338 |proj g|= 0.20161
## At iterate 20 f = -1338 |proj g|= 0.50961
## At iterate 30 f = -1338.1 |proj g|= 0.10069
## At iterate 40 f = -1338.1 |proj g|= 0.17518
## At iterate 50 f = -1338.1 |proj g|= 0.27559
## At iterate 60 f = -1338.1 |proj g|= 0.06003
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 69
## function evaluations 111
## segments explored during Cauchy searches 72
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00858394
## final function value -1338.13
##
## F = -1338.13
## final value -1338.131602
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1032.7 |proj g|= 20
## At iterate 10 f = -1335.3 |proj g|= 4.0886
## At iterate 20 f = -1338.1 |proj g|= 0.59767
## At iterate 30 f = -1338.1 |proj g|= 1.0666
## At iterate 40 f = -1338.1 |proj g|= 0.045435
## At iterate 50 f = -1338.1 |proj g|= 0.078825
##
## iterations 50
## function evaluations 66
## segments explored during Cauchy searches 55
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0788252
## final function value -1338.13
##
## F = -1338.13
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1338.131578
## converged
##
##
## ***************************
## Estimation of CV-set 3/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1322.9 |proj g|= 9.3055
## At iterate 10 f = -1323.3 |proj g|= 0.44233
## At iterate 20 f = -1323.4 |proj g|= 0.077529
##
## iterations 27
## function evaluations 46
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0242257
## final function value -1323.41
##
## F = -1323.41
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1323.412118
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1025.7 |proj g|= 20
## At iterate 10 f = -1319.3 |proj g|= 20.067
## At iterate 20 f = -1322.9 |proj g|= 0.17057
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 29
## function evaluations 74
## segments explored during Cauchy searches 35
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0414085
## final function value -1323.41
##
## F = -1323.41
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1323.412120
## converged
##
##
## ***************************
## Estimation of CV-set 4/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1321.3 |proj g|= 16.868
## At iterate 10 f = -1322.2 |proj g|= 0.14931
## At iterate 20 f = -1322.2 |proj g|= 0.48659
## At iterate 30 f = -1322.2 |proj g|= 0.019047
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1322.239329
## stopped after 32 iterations
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1019.9 |proj g|= 20
## At iterate 10 f = -1319.5 |proj g|= 3.5891
## At iterate 20 f = -1322.2 |proj g|= 0.033125
##
## iterations 22
## function evaluations 29
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0418117
## final function value -1322.23
##
## F = -1322.23
## final value -1322.233610
## converged
##
##
## ***************************
## Estimation of CV-set 5/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1346.7 |proj g|= 10.31
## At iterate 10 f = -1347.3 |proj g|= 0.24337
## At iterate 20 f = -1347.4 |proj g|= 0.78481
## At iterate 30 f = -1347.5 |proj g|= 0.48791
## At iterate 40 f = -1347.5 |proj g|= 0.076001
## At iterate 50 f = -1347.5 |proj g|= 0.053906
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 53
## function evaluations 93
## segments explored during Cauchy searches 54
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0640349
## final function value -1347.52
##
## F = -1347.52
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1347.517035
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1043.4 |proj g|= 20
## At iterate 10 f = -1345.7 |proj g|= 3.439
## At iterate 20 f = -1347.5 |proj g|= 0.2081
## At iterate 30 f = -1347.5 |proj g|= 1.2576
## At iterate 40 f = -1347.5 |proj g|= 0.13001
##
## iterations 46
## function evaluations 63
## segments explored during Cauchy searches 51
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00663204
## final function value -1347.52
##
## F = -1347.52
## final value -1347.517059
## converged
##
##
## ***************************
## Estimation of CV-set 6/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1320.9 |proj g|= 10.182
## At iterate 10 f = -1321.2 |proj g|= 0.27646
## At iterate 20 f = -1321.2 |proj g|= 0.066922
## At iterate 30 f = -1321.2 |proj g|= 0.035871
##
## iterations 30
## function evaluations 36
## segments explored during Cauchy searches 31
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0358709
## final function value -1321.18
##
## F = -1321.18
## final value -1321.183539
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1008.5 |proj g|= 20
## At iterate 10 f = -1319.2 |proj g|= 3.6121
## At iterate 20 f = -1321.2 |proj g|= 0.097296
##
## iterations 24
## function evaluations 30
## segments explored during Cauchy searches 29
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0223154
## final function value -1321.17
##
## F = -1321.17
## final value -1321.171547
## converged
##
##
## ***************************
## Estimation of CV-set 7/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1330.9 |proj g|= 2.1047
## At iterate 10 f = -1332 |proj g|= 1.8288
## At iterate 20 f = -1332.1 |proj g|= 0.66345
## At iterate 30 f = -1332.3 |proj g|= 0.14817
##
## iterations 38
## function evaluations 43
## segments explored during Cauchy searches 39
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.000450916
## final function value -1332.32
##
## F = -1332.32
## final value -1332.320536
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1017.7 |proj g|= 20
## At iterate 10 f = -1328.6 |proj g|= 2.9518
## At iterate 20 f = -1332.2 |proj g|= 0.43645
## At iterate 30 f = -1332.3 |proj g|= 0.45755
## At iterate 40 f = -1332.3 |proj g|= 0.021828
##
## iterations 43
## function evaluations 58
## segments explored during Cauchy searches 48
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0851882
## final function value -1332.32
##
## F = -1332.32
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1332.320410
## converged
##
##
## ***************************
## Estimation of CV-set 8/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1308.2 |proj g|= 12.34
## At iterate 10 f = -1309.5 |proj g|= 0.24255
## At iterate 20 f = -1309.6 |proj g|= 0.047146
## At iterate 30 f = -1309.6 |proj g|= 0.080695
## At iterate 40 f = -1309.6 |proj g|= 0.011327
##
## iterations 40
## function evaluations 58
## segments explored during Cauchy searches 40
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0113275
## final function value -1309.59
##
## F = -1309.59
## final value -1309.594791
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1003.1 |proj g|= 20
## At iterate 10 f = -1307.2 |proj g|= 2.8981
## At iterate 20 f = -1309.3 |proj g|= 1.2754
## At iterate 30 f = -1309.4 |proj g|= 1.6429
## At iterate 40 f = -1309.5 |proj g|= 0.22026
## At iterate 50 f = -1309.5 |proj g|= 1.8767
## At iterate 60 f = -1309.6 |proj g|= 0.094154
## At iterate 70 f = -1309.6 |proj g|= 0.012398
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 72
## function evaluations 122
## segments explored during Cauchy searches 78
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0509512
## final function value -1309.59
##
## F = -1309.59
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1309.594828
## converged
##
##
## ***************************
## Estimation of CV-set 9/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1353 |proj g|= 8.74
## At iterate 10 f = -1353.2 |proj g|= 0.21479
##
## iterations 13
## function evaluations 16
## segments explored during Cauchy searches 13
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0334872
## final function value -1353.24
##
## F = -1353.24
## final value -1353.238180
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1046.2 |proj g|= 20
## At iterate 10 f = -1350.5 |proj g|= 3.0881
## At iterate 20 f = -1353.2 |proj g|= 0.13529
##
## iterations 23
## function evaluations 30
## segments explored during Cauchy searches 28
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0461884
## final function value -1353.2
##
## F = -1353.2
## final value -1353.203269
## converged
##
##
## ***************************
## Estimation of CV-set 10/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1334.7 |proj g|= 9.2237
## At iterate 10 f = -1335.3 |proj g|= 0.31064
## At iterate 20 f = -1335.3 |proj g|= 0.1466
##
## iterations 25
## function evaluations 38
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0154309
## final function value -1335.29
##
## F = -1335.29
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1335.293048
## converged
##
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1021.3 |proj g|= 20
## At iterate 10 f = -1333.1 |proj g|= 2.4865
## At iterate 20 f = -1335.1 |proj g|= 0.69007
## At iterate 30 f = -1335.2 |proj g|= 0.98087
## At iterate 40 f = -1335.3 |proj g|= 0.024394
##
## iterations 40
## function evaluations 49
## segments explored during Cauchy searches 45
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0243941
## final function value -1335.29
##
## F = -1335.29
## final value -1335.292899
## converged
##
print(est.denver.B.cv)
## Cross-validation parameter estimation for STmodel
## with 10 CV-groups and 2 starting points.
## Results: 10 converged, 0 not converged.
##
## No fixed parameters.
##
## Estimated function values and convergence info:
## value convergence conv eigen.min eigen.all.min
## 1 1325.551 TRUE TRUE 0.3083541248 NA
## 2 1338.132 TRUE TRUE 0.0013147477 NA
## 3 1323.412 TRUE TRUE 1.0261104123 NA
## 4 1322.234 TRUE TRUE 0.0019628702 NA
## 5 1347.517 TRUE TRUE 0.0014229949 NA
## 6 1321.184 TRUE TRUE 0.0903619032 NA
## 7 1332.321 TRUE TRUE 0.0008020935 NA
## 8 1309.595 TRUE TRUE 0.1956246001 NA
## 9 1353.238 TRUE TRUE 0.2799606817 NA
## 10 1335.293 TRUE TRUE 0.5380999282 NA
par(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.B, pars="all"),
plotCI((1:length(par))+.3, par, uiw=1.96*sd,
col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.B.cv, "all", boxwex=.4, col="grey", add=TRUE)
#' Save the results as an .rdata object
save(denver.data.B, denver.model.B, est.denver.model.B, est.denver.B.cv,
file = here::here("Results", "Denver_ST_Model_B.rdata"))
Making predictions using the CV model. Printing out the CV summary statistics as well
pred.B.cv <- predictCV(denver.model.B, est.denver.B.cv, LTA = T)
pred.B.cv.log <- predictCV(denver.model.B, est.denver.B.cv, LTA = T,
transform="unbiased")
head(pred.B.cv$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 0.2256719 0.2715704 0.4040500 0.3325563 0.3381283 0.3552283
## 2009-01-05 NA 0.2454387 0.2867163 0.4125644 0.3424156 0.3470399 0.3569286
## 2009-01-12 NA 0.3202412 0.3522732 0.4671434 0.3973428 0.4010627 0.4036430
## 2009-01-19 NA 0.2611628 0.2956337 0.4101100 0.3424710 0.3453806 0.3404619
## 2009-01-26 NA 0.2625262 0.2945710 0.4039017 0.3370371 0.3392808 0.3265794
## 2009-02-02 NA 0.2774678 0.3063329 0.4094768 0.3426068 0.3443714 0.3235264
## d_15 d_16 d_17 d_18 d_19 d_2
## 2008-12-29 0.08474891 0.4246954 0.2028410 0.3308776 0.2007175 0.3439717
## 2009-01-05 0.12050759 0.4430645 0.2272220 0.3412991 0.2302263 0.3775853
## 2009-01-12 0.21196409 0.5071419 0.3034162 0.4055201 0.3213495 0.4671328
## 2009-01-19 0.16712856 0.4627590 0.2527525 0.3395599 0.2637587 0.4242432
## 2009-01-26 0.18303365 0.4698510 0.2585887 0.3330385 0.2743600 0.4444420
## 2009-02-02 0.21192172 0.4905654 0.2764793 0.3404257 0.3011394 0.4815097
## d_20 d_21 d_22 d_23 d_25 d_26
## 2008-12-29 0.3708507 0.2226343 0.2735033 0.3692873 0.2239682 0.2836085
## 2009-01-05 0.3803355 0.2426663 0.2926294 0.3848083 0.2417591 0.3000840
## 2009-01-12 0.4359096 0.3129022 0.3667768 0.4456929 0.3156517 0.3718684
## 2009-01-19 0.3799220 0.2604923 0.3065849 0.3973926 0.2541472 0.3099290
## 2009-01-26 0.3748357 0.2629673 0.3065892 0.3994628 0.2547692 0.3091289
## 2009-02-02 0.3816306 0.2773724 0.3197456 0.4137345 0.2701712 0.3226977
## d_28 d_29 d_3 d_31 d_32 d_33
## 2008-12-29 0.2037387 0.2646549 0.3040163 0.3402394 0.2963519 0.1975749
## 2009-01-05 0.2205081 0.2741179 0.3174851 0.3546011 0.3047909 0.2182837
## 2009-01-12 0.2933147 0.3338391 0.3759481 0.4245347 0.3669623 0.3005877
## 2009-01-19 0.2305899 0.2710521 0.3244478 0.3617353 0.2988150 0.2341560
## 2009-01-26 0.2297874 0.2633712 0.3221305 0.3611373 0.2898981 0.2358924
## 2009-02-02 0.2435023 0.2679122 0.3304864 0.3763881 0.2946226 0.2537815
## d_34 d_35 d_36 d_37 d_38 d_39
## 2008-12-29 0.3537773 0.3138491 0.3399919 0.3601579 0.2811222 0.3474737
## 2009-01-05 0.3678469 0.3254832 0.3488600 0.3679288 0.2988888 0.3520544
## 2009-01-12 0.4438244 0.3891428 0.4115710 0.4215859 0.3705696 0.4120839
## 2009-01-19 0.3717432 0.3264221 0.3441967 0.3632544 0.3122969 0.3391408
## 2009-01-26 0.3688713 0.3209440 0.3364089 0.3551737 0.3138032 0.3280427
## 2009-02-02 0.3834977 0.3284855 0.3427212 0.3581371 0.3296518 0.3323407
## d_4 d_40 d_41 d_42 d_43 d_44
## 2008-12-29 0.2382955 0.2492708 0.3870883 0.2337331 0.3137349 0.3119008
## 2009-01-05 0.2686225 0.2713930 0.3883383 0.2593849 0.3231533 0.3263064
## 2009-01-12 0.3808060 0.3754023 0.4356365 0.3668958 0.3865030 0.3960576
## 2009-01-19 0.2981158 0.2846203 0.3713036 0.2795488 0.3199534 0.3321672
## 2009-01-26 0.3073177 0.2858647 0.3577746 0.2841226 0.3132742 0.3295439
## 2009-02-02 0.3403762 0.3111443 0.3560063 0.3125931 0.3210603 0.3414555
## d_45 d_46 d_47 d_48 d_49 d_5
## 2008-12-29 0.3965320 0.2567529 0.1750927 0.3120921 0.3362977 0.2812456
## 2009-01-05 0.3963727 0.2742614 0.2026760 0.3250762 0.3459776 0.2983658
## 2009-01-12 0.4424060 0.3478701 0.2918032 0.3941431 0.4008662 0.3706735
## 2009-01-19 0.3771173 0.2860787 0.2320709 0.3277827 0.3462545 0.3094279
## 2009-01-26 0.3631061 0.2864096 0.2403103 0.3235076 0.3415375 0.3091272
## 2009-02-02 0.3614669 0.3015152 0.2644451 0.3339626 0.3484125 0.3230546
## d_50 d_51 d_52 d_53 d_54 d_55
## 2008-12-29 0.2118696 0.2859845 0.3952992 0.1619867 0.2032670 0.2547967
## 2009-01-05 0.2300339 0.2994959 0.4021042 0.1961947 0.2270731 0.2729522
## 2009-01-12 0.3043488 0.3650713 0.4549500 0.3121661 0.3125292 0.3728944
## 2009-01-19 0.2433689 0.3043440 0.3961372 0.2330581 0.2493588 0.2778413
## 2009-01-26 0.2446715 0.3009760 0.3880805 0.2455238 0.2545202 0.2745069
## 2009-02-02 0.2609549 0.3107765 0.3917193 0.2814335 0.2760431 0.2948136
## d_56 d_57 d_58 d_6 d_60 d_61
## 2008-12-29 0.2854202 0.4341711 0.2565010 0.2954464 0.2226554 0.2857551
## 2009-01-05 0.3015254 0.4326568 0.2746232 0.3090520 0.2443350 0.3015237
## 2009-01-12 0.3728553 0.4816223 0.3430821 0.3779759 0.3276838 0.3675802
## 2009-01-19 0.3107126 0.4085735 0.2891807 0.3132019 0.2624518 0.3111780
## 2009-01-26 0.3096394 0.3913963 0.2906024 0.3096097 0.2656235 0.3099508
## 2009-02-02 0.3229555 0.3874346 0.3045200 0.3204431 0.2852508 0.3210304
## d_62 d_63 d_64 d_65 d_66 d_67
## 2008-12-29 0.3313692 0.2867062 0.3873173 0.2832015 0.2454683 0.2485100
## 2009-01-05 0.3458785 0.3005263 0.3861333 0.3000520 0.2633762 0.2763734
## 2009-01-12 0.4156277 0.3663932 0.4308383 0.3707025 0.3364746 0.3862913
## 2009-01-19 0.3515107 0.3059206 0.3635800 0.3111593 0.2760240 0.3017601
## 2009-01-26 0.3483175 0.3027505 0.3466204 0.3110289 0.2765243 0.3097724
## 2009-02-02 0.3592168 0.3126757 0.3407722 0.3247696 0.2912599 0.3424822
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.2816600 0.2675658 0.2719639 NA NA NA NA NA
## 2009-01-05 0.2998819 0.2833328 0.2883306 NA NA NA NA NA
## 2009-01-12 0.3700401 0.3493903 0.3567725 NA NA NA NA NA
## 2009-01-19 0.3136151 0.2929949 0.2989293 NA NA NA NA NA
## 2009-01-26 0.3141158 0.2917834 0.2984695 NA NA NA NA NA
## 2009-02-02 0.3272234 0.3028901 0.3112074 NA NA NA NA NA
tail(pred.B.cv$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 1.0592974 1.0228608 1.061747 1.0395015 1.0352110 0.9430054
## 2020-11-09 NA 0.6743546 0.6436968 0.670476 0.6555533 0.6414596 0.5832082
## 2020-11-16 NA 1.0748031 1.0491068 1.086222 1.0686527 1.0601094 1.0035757
## 2020-11-23 NA 1.0128161 0.9815864 1.026519 1.0057695 0.9995560 0.9514267
## 2020-11-30 NA 1.2072519 1.1746121 1.236153 1.2102909 1.2097030 1.1687079
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2
## 2020-11-02 1.1419553 1.1275054 1.1310229 1.0238101 1.1039037 1.2837506
## 2020-11-09 0.7476813 0.7102194 0.7356171 0.6491306 0.7212576 0.8284186
## 2020-11-16 1.1288283 1.1099571 1.1315912 1.0667230 1.1094668 1.2016631
## 2020-11-23 1.0364911 1.0045859 1.0578360 1.0122618 1.0301734 1.0992655
## 2020-11-30 1.1975460 1.1997819 1.2452609 1.2195107 1.1992809 1.2551670
## 2020-12-07 NA NA NA NA NA NA
## d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 2020-11-02 1.0191267 1.0805214 1.1344515 1.0770611 0.9911015 1.063403 0.9644699
## 2020-11-09 0.6419187 0.6906482 0.7454757 0.6860287 0.6114598 0.662137 0.5896073
## 2020-11-16 1.0581193 1.0906546 1.1559272 1.0892098 1.0168345 1.078312 0.9941662
## 2020-11-23 0.9940127 1.0211557 1.0949061 1.0059062 0.9483883 1.018356 0.9300305
## 2020-11-30 1.2055783 1.2123088 1.2955459 1.1892691 1.1432666 1.222471 1.1271004
## 2020-12-07 NA NA NA NA NA NA NA
## d_29 d_3 d_31 d_32 d_33 d_34
## 2020-11-02 0.9756170 1.0777925 0.9988789 0.9793660 1.0054510 1.0214732
## 2020-11-09 0.6090801 0.6853581 0.6207785 0.6114122 0.6246790 0.6480789
## 2020-11-16 1.0240521 1.0964536 1.0265601 1.0334813 1.0267347 1.0591085
## 2020-11-23 0.9700332 1.0282825 0.9582950 0.9824494 0.9663557 1.0026145
## 2020-11-30 1.1776068 1.2287002 1.1473745 1.1940734 1.1596187 1.2034563
## 2020-12-07 NA NA NA NA NA NA
## d_35 d_36 d_37 d_38 d_39 d_4
## 2020-11-02 1.0415438 1.0048275 1.0566150 1.0367478 0.9030595 1.0685735
## 2020-11-09 0.6667445 0.6333597 0.6803382 0.6496136 0.5461058 0.6657496
## 2020-11-16 1.0822076 1.0529320 1.1005811 1.0523396 0.9632485 1.0695919
## 2020-11-23 1.0240669 0.9997173 1.0412693 0.9802670 0.9176003 1.0111316
## 2020-11-30 1.2313727 1.2084478 1.2506371 1.1743400 1.1272262 1.2000965
## 2020-12-07 NA NA NA NA NA NA
## d_40 d_41 d_42 d_43 d_44 d_45
## 2020-11-02 0.8413985 0.9707494 1.0133947 0.9667202 1.0469778 0.9281184
## 2020-11-09 0.4871106 0.6037283 0.6257857 0.5950959 0.6599757 0.5622223
## 2020-11-16 0.8729626 1.0334131 1.0379072 1.0126210 1.0738734 0.9901673
## 2020-11-23 0.8106980 0.9794114 0.9843721 0.9539198 1.0150174 0.9377863
## 2020-11-30 0.9991102 1.1997688 1.1778619 1.1613094 1.2177293 1.1555613
## 2020-12-07 NA NA NA NA NA NA
## d_46 d_47 d_48 d_49 d_5 d_50
## 2020-11-02 1.0114037 1.0703124 1.0200884 1.0144921 1.0151890 0.9595869
## 2020-11-09 0.6285872 0.6935142 0.6472952 0.6318189 0.6376615 0.5799605
## 2020-11-16 1.0296364 1.0842739 1.0581806 1.0464654 1.0321577 0.9765684
## 2020-11-23 0.9687930 1.0127037 1.0034845 0.9787747 0.9695431 0.9127071
## 2020-11-30 1.1631952 1.1797399 1.2050546 1.1762742 1.1639134 1.1030276
## 2020-12-07 NA NA NA NA NA NA
## d_51 d_52 d_53 d_54 d_55 d_56
## 2020-11-02 1.0180081 1.0472624 1.1533852 1.0267021 0.9534371 1.0170177
## 2020-11-09 0.6499887 0.6646837 0.7306014 0.6454562 0.5876489 0.6347523
## 2020-11-16 1.0606907 1.0842299 1.1325974 1.0453139 1.0116455 1.0350990
## 2020-11-23 0.9989803 1.0264275 1.0687857 0.9792002 0.9695331 0.9763999
## 2020-11-30 1.1899641 1.2363109 1.2482694 1.1734064 1.1756000 1.1742615
## 2020-12-07 NA NA NA NA NA NA
## d_57 d_58 d_6 d_60 d_61 d_62 d_63
## 2020-11-02 0.9611157 1.0753301 1.0284772 1.0142264 1.053816 1.0684598 1.0344341
## 2020-11-09 0.6116437 0.6854853 0.6472296 0.6319410 0.667495 0.6844248 0.6548452
## 2020-11-16 1.0404840 1.0871688 1.0638627 1.0287715 1.069091 1.0941749 1.0645456
## 2020-11-23 1.0018217 1.0186686 1.0041998 0.9629806 1.007879 1.0384102 1.0046810
## 2020-11-30 1.2233153 1.2082711 1.2071578 1.1537165 1.207289 1.2432936 1.2054039
## 2020-12-07 NA NA NA NA NA NA NA
## d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 0.9675467 1.0558200 1.0256708 1.0136834 1.1363678 1.0344538
## 2020-11-09 0.6091437 0.6490888 0.6429818 0.6305675 0.7437944 0.6484312
## 2020-11-16 1.0417461 1.0617888 1.0468392 1.0376018 1.1498478 1.0503399
## 2020-11-23 0.9985844 1.0027639 0.9801774 0.9728781 1.0854197 0.9874180
## 2020-11-30 1.2229546 1.2085839 1.1764166 1.1494601 1.2794957 1.1866676
## 2020-12-07 NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.0464963 NA NA NA NA NA
## 2020-11-09 0.6604396 NA NA NA NA NA
## 2020-11-16 1.0660222 NA NA NA NA NA
## 2020-11-23 1.0015201 NA NA NA NA NA
## 2020-11-30 1.1976984 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA
head(pred.B.cv.log$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 1.285715 1.346294 1.537529 1.432267 1.440270 1.465110
## 2009-01-05 NA 1.311077 1.366507 1.550255 1.446051 1.452754 1.467191
## 2009-01-12 NA 1.412620 1.458780 1.636827 1.527324 1.533016 1.536977
## 2009-01-19 NA 1.331348 1.378201 1.545771 1.445474 1.449686 1.442573
## 2009-01-26 NA 1.332971 1.376532 1.535948 1.437400 1.440628 1.422446
## 2009-02-02 NA 1.352884 1.392658 1.544336 1.445244 1.447797 1.417930
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2008-12-29 1.116118 1.570506 1.256995 1.428666 1.253541 1.448685 1.487322
## 2009-01-05 1.156607 1.599171 1.287688 1.443225 1.290773 1.497763 1.501089
## 2009-01-12 1.267234 1.704576 1.389325 1.538570 1.413621 1.637645 1.586493
## 2009-01-19 1.211556 1.630237 1.320438 1.440057 1.334269 1.568547 1.499804
## 2009-01-26 1.230881 1.641564 1.327962 1.430453 1.348290 1.600270 1.491947
## 2009-02-02 1.266877 1.675710 1.351773 1.440872 1.384725 1.660480 1.501925
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2008-12-29 1.281998 1.350112 1.485854 1.282840 1.363825 1.257149 1.337016
## 2009-01-05 1.307619 1.375774 1.508671 1.305704 1.386068 1.278249 1.349400
## 2009-01-12 1.402463 1.481275 1.602984 1.405684 1.488836 1.374634 1.432135
## 2009-01-19 1.330610 1.394438 1.527082 1.321707 1.399109 1.290935 1.344736
## 2009-01-26 1.333709 1.394198 1.529990 1.322424 1.397744 1.289796 1.334248
## 2009-02-02 1.352904 1.412474 1.551785 1.342863 1.416650 1.307524 1.340166
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2008-12-29 1.391968 1.441786 1.380182 1.249607 1.460869 1.404571 1.441747
## 2009-01-05 1.410446 1.462301 1.391486 1.275450 1.481214 1.420644 1.454179
## 2009-01-12 1.494994 1.567904 1.480375 1.384574 1.597802 1.513681 1.547908
## 2009-01-19 1.419655 1.472210 1.382561 1.295350 1.486416 1.421387 1.446750
## 2009-01-26 1.416132 1.471117 1.370055 1.297409 1.481935 1.413403 1.435282
## 2009-02-02 1.427833 1.493555 1.376364 1.320676 1.503596 1.423934 1.444183
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2008-12-29 1.471503 1.359322 1.452254 1.302283 1.316655 1.511670 1.296355
## 2009-01-05 1.482580 1.383297 1.458582 1.342038 1.345762 1.513150 1.329698
## 2009-01-12 1.563930 1.485725 1.548503 1.501024 1.492935 1.586060 1.480289
## 2009-01-19 1.475013 1.401327 1.439319 1.381632 1.363111 1.486934 1.356216
## 2009-01-26 1.462899 1.403201 1.423227 1.394185 1.364594 1.466708 1.362219
## 2009-02-02 1.467051 1.425431 1.429195 1.440870 1.399360 1.463928 1.401389
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2008-12-29 1.404384 1.402961 1.526013 1.325594 1.221827 1.401019 1.437636
## 2009-01-05 1.417273 1.422895 1.525356 1.348839 1.255698 1.419152 1.451211
## 2009-01-12 1.509587 1.525289 1.596833 1.451711 1.372464 1.520464 1.532715
## 2009-01-19 1.412098 1.430571 1.495604 1.364592 1.292652 1.422704 1.450953
## 2009-01-26 1.402459 1.426572 1.474549 1.364934 1.303154 1.416522 1.443883
## 2009-02-02 1.413237 1.443474 1.471944 1.385621 1.334834 1.431318 1.453660
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2008-12-29 1.359190 1.267413 1.365974 1.524133 1.206605 1.256740 1.323951
## 2009-01-05 1.382338 1.290483 1.384200 1.534124 1.248274 1.286710 1.347862
## 2009-01-12 1.485689 1.389885 1.477680 1.616990 1.401450 1.401207 1.489195
## 2009-01-19 1.397182 1.307538 1.390349 1.524322 1.294608 1.315193 1.353902
## 2009-01-26 1.396559 1.309137 1.385460 1.511839 1.310641 1.321804 1.349183
## 2009-02-02 1.415985 1.330544 1.398939 1.517154 1.358396 1.350406 1.376693
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2008-12-29 1.364876 1.584005 1.326158 1.380065 1.281344 1.365527 1.430542
## 2009-01-05 1.386712 1.581223 1.350082 1.398554 1.309114 1.386892 1.451018
## 2009-01-12 1.488934 1.660220 1.445434 1.497957 1.422603 1.481282 1.555433
## 2009-01-19 1.398978 1.542985 1.369336 1.403696 1.332526 1.399791 1.458513
## 2009-01-26 1.397275 1.516481 1.371080 1.398416 1.336562 1.397867 1.453607
## 2009-02-02 1.415845 1.510311 1.390135 1.413460 1.362897 1.413278 1.469341
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2008-12-29 1.366960 1.512016 1.362151 1.311422 1.315654 1.360079 1.340914
## 2009-01-05 1.385627 1.509817 1.384907 1.334807 1.352481 1.384735 1.361892
## 2009-01-12 1.479634 1.578468 1.485922 1.435739 1.509280 1.485040 1.454581
## 2009-01-19 1.392543 1.475494 1.399734 1.351281 1.386676 1.403299 1.374569
## 2009-01-26 1.387921 1.450439 1.399313 1.351762 1.397612 1.403785 1.372700
## 2009-02-02 1.401599 1.441796 1.418489 1.371673 1.443908 1.422138 1.387871
## d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 1.346955 NA NA NA NA NA
## 2009-01-05 1.368831 NA NA NA NA NA
## 2009-01-12 1.465467 NA NA NA NA NA
## 2009-01-19 1.382841 NA NA NA NA NA
## 2009-01-26 1.381992 NA NA NA NA NA
## 2009-02-02 1.399542 NA NA NA NA NA
tail(pred.B.cv.log$pred.all$EX)
## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 2.899503 2.794215 2.904608 2.842982 2.829757 2.581567
## 2020-11-09 NA 1.972851 1.912181 1.963776 1.936143 1.908337 1.801091
## 2020-11-16 NA 2.944390 2.867957 2.975923 2.926122 2.900149 2.741870
## 2020-11-23 NA 2.767585 2.680783 2.803544 2.747680 2.729643 2.602450
## 2020-11-30 NA 3.362062 3.251912 3.457825 3.371360 3.368125 3.234170
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 3.145770 3.101308 3.114198 2.796847 3.029272 3.627302 2.782459
## 2020-11-09 2.120700 2.042818 2.096767 1.922487 2.065863 2.300051 1.907834
## 2020-11-16 3.104689 3.046361 3.115183 2.918621 3.045665 3.340231 2.892463
## 2020-11-23 2.831069 2.741593 2.893688 2.763885 2.813628 3.014971 2.712944
## 2020-11-30 3.326239 3.332671 3.490455 3.400611 3.332469 3.523719 3.352551
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2020-11-02 2.960609 3.126655 2.949135 2.706198 2.910953 2.635964 2.666323
## 2020-11-09 2.004469 2.118615 1.994257 1.851257 1.948362 1.811850 1.847847
## 2020-11-16 2.990173 3.193375 2.984204 2.776691 2.953610 2.715368 2.798094
## 2020-11-23 2.789498 3.004181 2.745579 2.593196 2.781582 2.546880 2.651032
## 2020-11-30 3.377457 3.671754 3.298255 3.151575 3.411527 3.102075 3.262947
## 2020-12-07 NA NA NA NA NA NA NA
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2020-11-02 2.952940 2.727177 2.675206 2.746483 2.789847 2.846940 2.744197
## 2020-11-09 1.994034 1.868337 1.851283 1.876524 1.920267 1.956733 1.892365
## 2020-11-16 3.007579 2.803318 2.823133 2.805100 2.896363 2.964346 2.878585
## 2020-11-23 2.809275 2.618494 2.682643 2.640871 2.737406 2.796909 2.729371
## 2020-11-30 3.432821 3.163955 3.315126 3.204329 3.346730 3.441461 3.363126
## 2020-12-07 NA NA NA NA NA NA NA
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2020-11-02 2.889273 2.832845 2.480460 2.926076 2.330771 2.651326 2.768458
## 2020-11-09 1.982918 1.923129 1.735635 1.955587 1.635197 1.836534 1.878613
## 2020-11-16 3.018475 2.876513 2.633963 2.928462 2.405023 2.822165 2.836581
## 2020-11-23 2.844743 2.676456 2.516583 2.762262 2.259916 2.673893 2.688806
## 2020-11-30 3.507704 3.249937 3.103936 3.337194 2.728784 3.333471 3.263182
## 2020-12-07 NA NA NA NA NA NA NA
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2020-11-02 2.641400 2.864568 2.541062 2.763561 2.929324 2.786918 2.770918
## 2020-11-09 1.821192 1.944860 1.762137 1.884506 2.009418 1.919577 1.889473
## 2020-11-16 2.764654 2.941595 2.703134 2.814360 2.970015 2.895074 2.860008
## 2020-11-23 2.607003 2.773316 2.565269 2.648435 2.765014 2.741189 2.672719
## 2020-11-30 3.208040 3.396618 3.189805 3.217173 3.268114 3.353800 3.256439
## 2020-12-07 NA NA NA NA NA NA NA
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 2.774650 2.624441 2.780413 2.862679 3.182716 2.804348 2.607647
## 2020-11-09 1.901951 1.795357 1.924008 1.952325 2.085071 1.915152 1.808526
## 2020-11-16 2.821729 2.669341 2.900926 2.969835 3.116602 2.856557 2.763377
## 2020-11-23 2.650626 2.504396 2.727318 2.803130 2.924035 2.673937 2.649507
## 2020-11-30 3.219768 3.029810 3.301504 3.458176 3.499305 3.247513 3.256185
## 2020-12-07 NA NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 2.779868 2.628354 2.944689 2.811351 2.769972 2.883928 2.928094
## 2020-11-09 1.896521 1.852884 1.993747 1.919744 1.889710 1.959500 1.993898
## 2020-11-16 2.830183 2.844904 2.979170 2.911560 2.810089 2.927741 3.003283
## 2020-11-23 2.668997 2.737092 2.782011 2.742786 2.631290 2.753983 2.840241
## 2020-11-30 3.253423 3.416094 3.363174 3.360054 3.184646 3.362102 3.486146
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 2.827670 2.643929 2.887731 2.801725 2.767600 3.130032 2.827912
## 2020-11-09 1.934203 1.847262 1.922346 1.910622 1.886487 2.113405 1.922014
## 2020-11-16 2.913377 2.846945 2.904164 2.861252 2.834017 3.171711 2.872630
## 2020-11-23 2.744084 2.726770 2.737672 2.676897 2.656486 2.973804 2.697526
## 2020-11-30 3.354308 3.413058 3.363550 3.257765 3.169908 3.611032 3.292651
## 2020-12-07 NA NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.861798 NA NA NA NA NA
## 2020-11-09 1.944927 NA NA NA NA NA
## 2020-11-16 2.917491 NA NA NA NA NA
## 2020-11-23 2.735245 NA NA NA NA NA
## 2020-11-30 3.328343 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA
summary(pred.B.cv)
## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 712 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX
## obs 0.13672884 0.13672884 0.10865488
## average 0.08203204 0.08203204 0.06010949
##
## R2:
## EX.mu EX.mu.beta EX
## obs 0.7318650 0.7318650 0.8306709
## average 0.7028679 0.7028679 0.8404603
##
## Coverage of 95% prediction intervals:
## EX
## obs 0.9550562
## average 0.9206349
summary(pred.B.cv.log)
## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 712 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.19696507 0.19696507 0.15055972 0.15053352
## average 0.09884076 0.09884076 0.06180058 0.06231628
##
## R2:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.6975357 0.6975357 0.8232688 0.8233303
## average 0.7160000 0.7160000 0.8889723 0.8871116
##
## Coverage of 95% prediction intervals:
## EX.pred
## obs 0.9550562
## average 0.9206349
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModB.jpeg"),
width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()
## quartz_off_screen
## 2
What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites
x.B <- coef(est.denver.model.B, pars = "cov")$par
x.B
## [1] -8.693332 -7.074670 -7.568345 13.812075 -3.008518 -4.818215
E.B <- predict(denver.model.B, est.denver.model.B, LTA = T,
transform="unbiased", pred.var=FALSE)
print(E.B)
## Prediction for STmodel.
##
## Regression parameters:
## 0 Spatio-temporal covariate(s).
## 24 beta-fields regression parameters in x$pars.
##
## Regression parameters are assumed to be known and
## prediction variances do NOT include
## uncertainties in regression parameters.
##
## Prediction of beta-fields, (x$beta):
## List of 2
## $ mu: num [1:69, 1:3] 0.339 0.182 0.183 0.225 0.19 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX: num [1:69, 1:3] 0.344 0.185 0.181 0.215 0.187 ...
## ..- attr(*, "dimnames")=List of 2
##
## Predictions for log-Gaussian field of type: unbiased
##
## Predictions for 624 times at 69 locations.
## List of 5
## $ EX.mu : num [1:624, 1:69] 1.45 1.48 1.6 1.53 1.55 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.mu.beta: num [1:624, 1:69] 1.49 1.52 1.64 1.56 1.58 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX : num [1:624, 1:69] 1.52 1.56 1.68 1.6 1.62 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.pred : num [1:624, 1:69] 1.53 1.56 1.69 1.61 1.62 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.EX : num [1:624, 1:69] 0.396 0.419 0.493 0.445 0.456 ...
## ..- attr(*, "dimnames")=List of 2
##
## Variances have been computed.
## List of 2
## $ log.VX : num [1:624, 1:69] 0.05 0.0499 0.0498 0.0497 0.0496 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.VX.pred: num [1:624, 1:69] 0.0581 0.058 0.0579 0.0578 0.0577 ...
## ..- attr(*, "dimnames")=List of 2
##
## Mean squared prediciton errors NOT been computed.
##
## 69 temporal averages have been compute.
## List of 1
## $ LTA:'data.frame': 69 obs. of 4 variables:
week_preds <- as.data.frame(E.B$EX) %>%
mutate(week = as.Date(rownames(E.B$EX)),
year = year(as.Date(rownames(E.B$EX))))
week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]
box_preds <- week_preds %>%
select(week, all_of(week_sites)) %>%
#filter(week %in% week_weeks) %>%
mutate(month = month(week),
year = year(week)) %>%
filter(year %in% c(2009:2020))
box_data <- box_preds %>%
pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)
## week month year location
## Min. :2009-01-05 Min. : 1.000 Min. :2009 Length:42364
## 1st Qu.:2011-12-26 1st Qu.: 4.000 1st Qu.:2011 Class :character
## Median :2014-12-22 Median : 7.000 Median :2014 Mode :character
## Mean :2014-12-22 Mean : 6.506 Mean :2014
## 3rd Qu.:2017-12-18 3rd Qu.: 9.000 3rd Qu.:2017
## Max. :2020-12-07 Max. :12.000 Max. :2020
##
## pred
## Min. :0.3674
## 1st Qu.:1.1699
## Median :1.4733
## Mean :1.5901
## 3rd Qu.:2.0054
## Max. :4.0348
## NA's :1080
#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
group_by(month) %>%
summarize(median = median(pred, na.rm=T)) %>%
arrange(desc(median))
box_summary
pred_box_plot <- ggplot(box_data) +
geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
show.legend = F) +
#scale_color_viridis(name = "Month", discrete = T) +
facet_wrap(. ~ year, scales = "free") +
xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
# theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot
| model | trend_data | cov_beta0 | cov_beta1 | cov_beta2 | cov_nu_error | no_basis_fns | df_per_year | all_converge | cv_rmse_obs | cv_rmse_avg | cv_r2_obs | cv_r2_avg | cv_coverage_obs | cv_coverage_avg |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Model A | BC + NO2 + PM2.5 | iid | iid | NA | exp | 1 | 4 | TRUE | 0.16 | 0.06 | 0.81 | 0.91 | 0.95 | 0.95 |
| Model B | BC + NO2 + PM2.5 | iid | iid | iid | exp | 2 | 4 | TRUE | 0.15 | 0.06 | 0.82 | 0.89 | 0.96 | 0.92 |